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Game control methods in avoidance 
of ships collisions 


Józef Lisowski, Prof. 
Gdynia Maritime University 


ABSTRACT 


The paper introduces application of selected methods of a game theory for automation of 
the processes of moving marine objects, the game control processes in marine navigation 
and the base mathematical model of game ship control. State equations, control and state 
constraints have been defined first and then control goal function in the form of payments 
— the integral payment and the final one. Multi-stage positional and multi-step matrix, 
non-cooperative and cooperative, game and optimal control algorithms in a collision 
situation has been presented. The considerations have been illustrated as an examples of 


a computer simulations mspg.12 and msmg.12 algorithms to determine a safe own ship 5 trajectory in the 
process of passing ships encountered in Kattegat Strait. 


Keywords: marine navigation; safety at sea; safe ship control; optimal control; differential games; 
positional games; matrix games; computer simulation 


INTRODUCTION 


The control of the ship’s movement may be treated as 
a multilevel problem shown in Fig. 1, which results from the 
division of the entire control system of ship - within the frame 
of the performance of the cargo carriage by the ship’s operator 
- into clearly determined subsystems which are ascribed 
appropriate layers of control. 


Route 
Planning 
Safe 
navigation 


Adaptive 
control 


: 
Y 


Fig. 1. Multilevel ship movement control system 


This is connected both with a large number of dimensions of 
the control vector and of the status of the process, its random, 
fuzzy and decision making characteristics - which are affected 
by strong interference generated by the current, wind and the 
sea wave motion on the one hand, and a complex nature of 
the equations describing the ship’s dynamics with non-linear 
and non-stationary characteristics. The determination of the 
global control of the steering systems has in practice become 
too costly and ineffective. 

The integral part of the entire system is the process of 
the ship’s movement control, which may be described with 
appropriate differential equations of the kinematics and 
dynamics ofa ship being an object of the control under a variety 
of the ship’s operational conditions such as: 

- stabilisation of the course or trajectory, 

- adjustment of the ship’s speed, 

- precise steering at small speeds in port with thrusters or 
adjustable-pitch propeller, 

- stabilisation of the ship’s rolling, 

- commanding the towing group, 

- dynamic positioning of the ship. 


The functional draft of the system corresponds to a certain 
actual arrangement of the equipment. The increasing demands 
with regard to the safety of navigation are forcing the ship’s 
operators to install the systems of integrated navigation on 
board their ships. By improving the ship’s control these systems 
increase the safety of navigation of a ship - which is a very 
expensive object of the value, including the cargo, and the 
effectiveness of the carriage goods by sea [1, 2]. 
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SAFE SHIP CONTROL IN COLLISION 
SITUATIONS 


The challenge in research for effective methods to prevent 
ship collisions has become important with the increasing size, 
speed and number of ships participating in sea carriage. An 
obvious contribution in increasing safety of shipping has been 
firstly the application of radars and then the development of 
ARPA (Automatic Radar Plotting Aids) anti-collision system 
[3, 4, 5]. 

The ARPA system enables to track automatically at least 20 
encountered j objects as is shown in Fig. 2, determination of 
their movement parameters (speed Vj, course Vj) and elements 
of approach to the own ship (Dhin = DCPA, - Distance of the 
Closest Point of Approach, T = TCPA ,- Time to the Closest 
Point of Approach) and also the assessment of the collision 
risk r; [6, 7, 8]. 


Fig. 2. Navigational situation representing the passing 
of the own ship with the j-th ship 


The risk value (1) is possible to define by referring the 
current situation of approach, described by parameters 
D} and T} p to the assumed evaluation of the situation as 
sate, determined by a safe distance of approach D, and a safe 
time T, — which are necessary to execute a collision avoiding 
manoeuvre with consideration of distance D; to j-th met ships 
(Fig. 3) [9, 10, 11]. 


Fig. 3. The ships collision risk space in a function 
of relative distance and time of approaching the j-th ships 
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1 
, 2 = 42 2| 2 
J. T. D. 
r= |wi] 2m | 4w, Tmin (Ps) (1) 
D, T, 


The weight coefficients w, and w, are depended on the state 
visibility at sea, dynamic length Ly and dynamic beam B, of 
the ship, kind of water region and in practice are equal: 


0<[w(L, Ba), Wa(La Ba)] < 1 2) 
La = 1.1(1 + 0.345V16) (3) 
B, = 1.1 + 0.767V) (4) 


The functional scope ofa standard ARPA system ends with 
the trial manoeuvre by altering the course +Ay or the ship’s 
speed +AVselected by the navigator as is shown in Fig. 4. 


HDGO95.4 "GYRO 
spp+0.20 kt WI LOGI 
TRACKPILOT FAILURE GHU 
SET 0.0 * 0.50 NW 


T N €OGI29.0 ' 
Le ILON 019:33.527 E 5000.00 ki 


cse | n7 spo izm 


0.03 W 
BcT si 


CLOCK: 00:00:00 
Apr 13 09:24:48 2007 ZT 


Fig. 4. The screen of SAM Electronics ARPA 
on the sailing vessel s/v DAR MLODZIEZY 


The problem of selecting such a manoeuvre is very difficult 
as the process of control is very complex since it is dynamic, 
non-linear, multi-dimensional, non-stationary and game making 
in its nature. 

In practice, methods of selecting a manoeuvre assume 
a form of appropriate steering algorithms supporting navigator 
decision in a collision situation. Algorithms are programmed 
into the memory of a Programmable Logic Controller PLC. 
This generates an option within the ARPA anti-collision system 
or a training simulator [12, 13, 14]. 

There are various methods for the avoidance of ships 
collision. The simplest method is determination of the 
manoeuvre of a change of course or a speed of own ship in 
relation to the most dangerous ship encountered. A more 
effective method is to determine safe trajectory of the ship 
[15, 16, 17, 18, 19, 20]. Most adequate to the real character 
of control process is determination of a game trajectory of the 
ship [21, 22, 23]. 


GAME CONTROL IN MARINE 
NAVIGATION 


The classical issues of the theory of the decision process 
in marine navigation include the safe steering of a ship. The 
problem of non-collision strategies in the steering at sea 
appeared in the Isaacs’ work [24] called “the father of the 
differential games” and was developed by many authors both 
within the context of the game theory [25, 26, 27, 28, 29] and 
also in the steering under uncertainty conditions [30]. 
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The definition of the problem of avoiding a collision 
seems to be quite obvious, however, apart from the issue 
of the uncertainty of information which may be a result of 
external factors (weather conditions, sea state), incomplete 
knowledge about other objects and imprecise nature of the 
recommendations concerning the right of way contained in 
International Regulations for Preventing Collision at Sea 
COLREG [6]. 

The problem of determining safe strategies is still an urgent 
issue as aresult of an ever increasing traffic of vessels on 
particular water areas. It is also important due to the increasing 
requirements as to the safety of shipping and environmental 
protection, from one side, and to the improving opportunities 
to use computer supporting the navigator’s duties. In order to 
ensure safe navigation the ships are obliged to observe legal 
requirements contained in the COLREG Rules. 

However, these Rules refer exclusively to two ships under 
good visibility conditions, in case of restricted visibility the 
Rules provide only recommendations of general nature and 
they are unable to consider all necessary conditions of the 
real process. Therefore the real process of the ships passing 
exercises occurs under the conditions of indefiniteness and 
conflict accompanied by an imprecise co-operation among the 
ships in the light of the legal regulations. 

Consequently, it is reasonable for ship operational purposes, 
to present this process and to develop and examine methods 
for a safe steering of the ship by applying the rules of the 
game theory. 

Anecessity to consider simultaneously the strategies of the 
encountered objects and the dynamic properties of the ships as 
the steering objects is a good reason for the application of the 
differential game model, often called the dynamic game, for 
the description of the processes [31]. 


PROCESSES OF GAME SHIP CONTROL 


Assuming that the dynamic movement of the ships in time 
occurs under the influence of the appropriate sets of steering: 


Dhin = min D,(t) 2 D, (5) 
where 
se a set of the own ship’s strategies, 
se) — aset of the j-th ship’s strategies, 
é€=0 -— denotes course and trajectory stabilisation, 
g= — denotes the execution of the anti-collision 
manoeuvre in order to minimize the risk of collision, 
which in practice is achieved by satisfying the 
following inequality: 
€ € 
Ee gl © 
Diy — the smallest distance of approach of the own ship 
and the j-th encountered object, 
D, — safe approach distance in the prevailing conditions 
depends on the visibility conditions at sea, the 
COLREG Rules and the ship’s dynamics. 
D, — current distance to the j-th object taken from the ARPA 
anti-collision system. 
e=- 1l -— refers to the manoeuvring of the ship in order to 


achieve the closest point of approach, for example 
during the approach of a rescue vessel, transfer of 
cargo from ship to ship, destruction the enemy’s 
ship, etc.). 


In the adopted describing symbols we can discriminate the 
following type of steering ship in order to achieve a determined 
goal: 


- basic type of steering, stabilization of the course or 
trajectory: [S®, s®] 


- avoidance of a collision by executing: 
a) own ship’s manoeuvres: [S(, s®] 
b) manoeuvres of the j-th ship: [S®, s®] 
c) co-operative manoeuvres: [S®, s] 
- encounter of the ships: Ise”, sc 
- situations of a unilateral dynamic game: (sO, s®)] and 
[SaF 


Dangerous situations resulting from a faulty assessment 
of the approaching process by one of the party with the 
other party’s failure to conduct observation - one ship is 
equipped with a radar or an anti-collision system, the other 
with a damaged radar or without this device [32]. 

- chasing situations which refer to a typical conflicting 
dynamic game: (sc? s] and [s@, Se . 


The first case usually represents regular optimal control, 
the second and third are unilateral games while the fourth and 
fifth cases represent the conflicting games. 


BASE MATHEMATICAL MODEL OF GAME 
SHIP CONTROL 


As the process of steering the ship in collision situations, 
when a greater number of objects is encountered, often occurs 
under the conditions of indefiniteness and conflict, accompanied 
by an inaccurate co-operation of the ships within the context 
of COLREG Regulations then the most adequate model of the 
process which has been adopted is a model of a dynamic game, 
in general of j tracked ships as objects of steering. 

The diversity of selection of possible models directly 
affects the synthesis of the ship’s handling algorithms which 
are afterwards effected by the ship’s handling device directly 
linked to the ARPA system and, consequently, determines the 
effects of the safe and optimal control. 

The most general description of the own ship passing the 
j number of other encountered ships is the model of a differential 
game of aj number of objects, shown in Fig. 5. 


Mic,  Uji"icj UmiU mem 


Lk KE KK 


1 j 


X41 ma | | "ma [rdn 
| Xi xi dj 


ul A — 0,1 
o1 own ship 

x 
"ot O h 


Fig. 5. Block diagram of a base dynamic game model 


m ship 


The properties of the process are described by the state 
equation: 


i=1l,..., Gd; +d.);j=1,....m (7) 


where: 

Xog, (t) — d, dimensional vector of the process state of the own 
_ ship determined in a time span t € [tọ tą], 

Xj; (t) — d, dimensional vector of the process state for the j-th 


ship, 
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Ug.c,(t)— c, dimensional control vector of the own ship, 
u ic) —G dimensional control vector of the j-th ship. 


Taking into consideration the equations reflecting the 
own ship’s hydromechanics and equations of the own ship’s 
movement relative to the j-th encountered ship, the equations 
of the general state of the process (7) take the form (8). 


Xo = Xo2 


Xo 2 = AyXg.2Xo3 $A2XKo 3 [Xo 3] + O1X_ 3/Xo.3/Uo 


X03 = 84X93 /Xo,3]Xo,4|Xo,40 +Xo4)+ 


TASK 4 2X0,3X0,4 eee FA7Xo3/Xo 3/7 


+4gX 4 5[Xo [Xo 6 +O2% 6 3X0 4/Xo3U 1 


Xo 4 = 43% q 3X0 4 A4X 6 3X04 


Ko 4|t8sXo2Xo2 T 
+a9X_ > +b 2X, 3U5) ®) 
Xo,5 = 10X05 + b3u, 2 
Xo 6 =A Xo6 tbyu,, 


xX: 


il = X43 +Xj 2X09 +X ji 3 COSX i 3 


Kio =—Xo2Xj, +X), 8inx,, 
X 3 =—Xo2 +4, )X,34j1 
Xj4 Samid +b; jujo 


The state variables are represented by the following 

values: 

X, = Y — course of the own ship, 

a a angular turning speed of the own ship, 

X,3 = V — speed of the own ship, 

X,4= — drift angle of the own ship, 

X,5 =n — rotational speed of the screw propeller of the own 
ship, 

X= H — pitch of the adjustable propeller of the own ship, 

x1 = Dj- distance to j-th object, or x, — its coordinate, 

Xio= N; — bearing of the j-th object, or y, — its coordinate, 

Xj 3 = Yj — course of the j-th object, or B — relative meeting 
angle, 

Xj 4 = Vj — speed of the j-th object, 

where: d, = 6, d; = 4. 


While the control values are represented by: 

u, ; = &, — reference rudder angle of the own ship, or \- angular 
turning speed of the own ship, or y - course of the 
own ship, depending of a kind approximated model 
of process, 

u, 2 =n, — reference rotational speed of the own ship’s screw 
propeller, or force of the propeller thrust of the own 
ship, or speed of the own ship, 

u,3 = H,- reference pitch of the adjustable propeller of the 

own ship, 

= yj — course of the j-th object, or wy; - angular turning speed 

of the j-th object, 

uj =V- speed of the j-th object, or force of the propeller 

: thrust of the j-th object, 

where: c, = 3, ¢; = 2. 


Uii 
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Values of coefficients of the process state equations (8) for 
the 12 000 DWT container ship are given in Table 1. 


Tab. 1. Coefficients of basic game model equations. 


Coefficient 


Measure 
- 4.143-10- 
1.858-10 
- 6.934:107 
- 3.177-10? 
- 4.435 
- 0.895 
- 9.284-10 
1.357-10 
0.624 
- 0.200 
- 0.100 
- 7.979:104 
1.134-10? 
- 1.554-10 
0.200 
0.100 
- 3.333-10% 
9.536-10? 


In example for j = 20 objects the base game model is 
represented by i = 86 state variables of process control. 

The constraints of the control and the state of the process 
are connected with the basic condition for the safe passing of 
the objects at a safe distance D, in compliance with COLREG 
Rules, generally in the following form: 


The constraints referred to as the ships domains in the 
marine navigation, may assume a shape of a circle, ellipse, 
hexagon, or parabola and may be generated for example by an 
artificial neural network as is shown in Fig. 6 [33, 34]. 

The synthesis of the decision making pattern of the object 
control leads to the determination of the optimal strategies of 
the players who determine the most favourable, under given 
conditions, conduct of the process. For the class of non- 
coalition games, often used in the control techniques, the most 
beneficial conduct of the own control object as a player with 
j-th object is the minimization of her goal function in the form 
of the payments — the integral payment and the final one: 


tk 


= fixa Ol dt+ r(t) +d(t,) > min (10) 


to 

The integral payment represents loss of way by the ship 
while passing the encountered objects and the final payment 
determines the final risk of collision r(t) relative to the j- 
th object and the final deflection of the ship d(t,) from the 
reference trajectory. 

Generally two types of the steering goals are taken into 
consideration - programmed steering u,(t) and positional 
steering u,[x,(t)]. The basis for the decision making steering 
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are the decision making patterns of the positional steering 
processes, the patterns with the feedback arrangement 
representing the dynamic games. 

The application of reductions in the description of the own 
ship’s dynamics and the dynamic of the j-th encountered ship 
and their movement kinematics lead to synthesis of game ship 
control algorithms in collisions situations. 


Fig. 6. The shapes of the neural domains in the situation 
of 20 encountered ships in Kattegat Strait 


ALGORITHMS OF GAME SHIP CONTROL 


Multi-stage positional game control algorithms 


The general model of dynamic game is simplified to the 
multi-stage positional game of j participants not co-operating 
among them [35, 36]. 

State variables and control values are represented by: 


Xo1 = Xo, Xo2 = Yo, Xj =X; X25 Yj 


us =V, Ugo = V, uji = Y; Uj — V; (11) 


The essence of the positional game is to subordinate the 
strategies of the own ship to the current positions p(t,) of the 
encountered objects at the current step k. In this way the process 
model takes into consideration any possible alterations of the 
course and speed of the encountered objects while steering is 
in progress. The current state of the process is determined by 
the co-ordinates of the own ship’s position and the positions 
of the encountered objects: 


X, = (Xo Yo); x= (X;, Y)J =1,2,....m (12) 


The system generates its steering at the moment t, on the 
basis of data received from the ARPA anti-collision system 
pertaining to the positions of the encountered objects: 


Xo(tk) 
X(t) 


It is assumed, according to the general concept of a multi- 
stage positional game, that at each discrete moment of time t, 
the own ship knows the positions of the objects. 


p(t, )= j=l, 2y: M > k=1, 2i K (13) 


The constraints for the state co-ordinates: 


{x,(0, OEP (14) 
are navigational constraints, while steering constraints: 
u, E So uj € S;J= 1,2, as m (15) 


take into consideration: the ships’ movement kinematics, 
recommendations of the COLREG Rules and the condition to 
maintain a safe passing distance as per relationship (6). 

The closed sets Soj and Soy? defined as the sets of 
acceptable strategies of the participants to the game towards 
one another: 


{S, [pO], S LOL (16) 


are dependent, which means that the choice of steering u, by 
the j-th object changes the sets of acceptable strategies of 
other ships. 


Multi-stage non-cooperative positional game control 
algorithm mspg_nc.12 


The optimal steering of the own ship u(t), equivalented for 
the current position p(t) to the optimal positional control už (p). 
The sets of acceptable strategies U. ,[p(t,)] are determined for 
the encountered ships relative to the own ship and initial sets 
U, [p(t,)] of acceptable strategies of the own ship relative 
to each one of the encountered ship. The pair of vectors uj 
and u,,, relative to each j-th ship is determined and then the 
optimal positional strategy for the own ship uğ (p) from the 
condition (10). 


ty 
I= min max min | uO dt= Sx, L) mer 


m al 7 
z UES, u ES «ty 
uE NS j= Sw Yaj 
j=l 


The function S, refers to the continuous function of the 
manoeuvring goal of the own ship, characterising the distance 
of the ship at the initial moment t, to the nearest turning point 
L, on the reference p,(t,) route of the voyage. 

The optimal control of the own ship is calculated at each 
discrete stage of the ship’s movement by applying the Simplex 
method to solve the problem of the triple linear programming, 
assuming the relationship (17) as the goal function and the 
control constraints (9). 


Multi-stage cooperative positional game control 
algorithm mspg_c.12 


The quality index of control for a cooperative game has 
the form: 


t 
I= min min min | u(t) dt= Sh (x, Ly) 


m N] 
, ueS u ES .to 
WE NS yj i ojn Yoj 
H 


(18) 


Multi-step matrix game control algorithms 


When leaving aside the ship’s dynamics equations the 
general model of a dynamic game for the process of preventing 
collisions is reduced to the matrix game of j participants non- 
co-operating among them [9]. 

The state and steering variables are represented by the 
following values: 


Xj 1 = Dj, X2 = N} Uy) = Y, U2 = V, Uy 1 = Yj Y2 = Vj 


j=1,2,...,.m oa. 
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The game matrix R = [r,(u;, u,)] includes the values of the 
collision risk r, determined from relation (1) on the basis of data 
obtained from the ARPA anti-collision system for the acceptable 
strategies u, of the own ship and acceptable strategies uj ofany 
particular number of j encountered objects. 

The problem of determining an optimal strategy may be 
reduced to the task of solving dual linear programming Simplex 
method. Mixed strategy components express the distribution 
of probability p,(u,, u;) of using pure strategies by the players 
(20, 21]. 


Multi-step non-cooperative matrix game control 
algorithm msmg_nc.12 


As a result of using the following form for the control goal: 


* . 

I, = min max r; 

i 
Uo Uj 


(20) 


the probability matrix P = [p; (u,,u,)] of using particular pure 
strategies may be obtained. 

The solution for the control problem is the strategy 
representing the highest probability: 


Ug = Ug{[P; (4,50; max} 


Multi-step cooperative matrix game control 
algorithm msmg_c.12 


(21) 


The quality index of control for a cooperative game has 
the form: 


* . . 
I, = min min r} 


Uo uj 


(22) 


R ai 


Fig. 7. The place of identification of navigational situation in Kattegat Strait 
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Fig. 8. The 12 minute speed vectors of own ship 
and 15 encountered ships in navigational situation in Kattegat Strait 
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Fig. 9. Computer simulation of mspg_ne.12 algorithm for safe manoeuvring 
of the own ship in situation of passing 15 encountered ships, D,=1.0 nm, 
d(t,)=2.49 nm (nautical mile) 
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Fig. 10. Computer simulation of mspg_c.12 algorithm for safe manoeuvring 
of the own ship in situation of passing 15 encountered ships, D,=1.0 nm, 
A(t,)=1.78 nm 
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Fig. 11. Computer simulation of msmg_ne.12 algorithm for safe 
manoeuvring of the own ship in situation of passing 15 encountered ships, 
D.=1.0 nm, d(t,)=1.86 nm 
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Fig. 12. Computer simulation of msmg_c.12 algorithm for safe 
manoeuvring of the own ship in situation of passing 15 encountered ships, 
D=1.0 nm, d(t,)=1.24 nm 


COMPUTER SIMULATION 


Computer simulation of control game algorithms were 
carried out on an example of a real navigational situations of 
passing j=15 encountered ships. The situations were registered 
in Kattegat Strait on board r/v HORYZONT II, a research 
and training vessel of the Gdynia Maritime University, on 
the radar screen of the ARPA anti-collision system Raytheon 
(Fig. 7 and 8). 

Examples of safe positional game trajectories are shown 
in Fig. 9 and 10. 

Examples of safe matrix game trajectories are shown in 
Fig. 11 and 12. 


CONCLUSIONS 


- The application of the models of a game theory for the 
synthesis of an optimal manoeuvring makes it possible 
to determine the safe game trajectory of the own ship 
in situations when she passes a greater number of the 
encountered ships. 

- Developed algorithms takes also into consideration the 
Rules of the COLREGS Rules and the advance time of 
the manoeuvre approximating ship’s dynamic properties 
and evaluates the final deviation of the real trajectory from 
reference value. 

- The positional game control algorithms determine game 
and safe trajectory of the own ship with relation to of all 
encountered ships. 

- The matrix game control algorithms determine game and 
safe trajectory of the own ship with relation to of the ship 
of most dangerous. 

- To sum up it may be stated that the control methods 
considered in this study are, in a certain sense, formal 
models for the thinking processes of a navigating 
officer steering of own ship and making decisions on 
manoeuvres. 

- Therefore they may be applied in the construction ofa new 
model of ARPA system containing a computer supporting 
the navigator’s decision making. 
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ABSTRACT 


The paper presents the physical model of the sea-going ship. The wooden 
hull of this ship was obtained from Ship Design and Research Centre 
and next rebuilt and reconstructed for future tests of the different control 
systems in the real-time experiments. The paper describes the succeeding 
stages of the work, steering and driving devices installed on board, power 
electric diagrams and navigational equipment designated for ship. The 
range of the possible applications of the constructed vessel are presented 


at the end of the paper 


Keywords: ship’s hull; main engine; tunnel thrusters; electrical motors; 
marine navigation devices; real-time experiments 


INTRODUCTION 


The seagoing cargo ship can be considered as a complex 
control object with its properties very different from 
manoeuvring inland vehicles or aircraft. One can distinguish, 
among others, three systems: 

- the navigation one with all measurement devices e.g., 
gyrocompass, log, echo sounder, radar, ECDIS. AIS, etc, 

- the ship’s movement steering one with subsystems e.g. 
autopilot, DSP, speed and/or trajectory controller etc, 

- and the power one with energy generation and distribution 
processes. 


The dissimilarity of the ship motion process in comparison 
to the adequate ones for inland vehicles or aircrafts is related 
to the following aspects: 

- very small power of the driving devices. The coefficient 
power/mass is significantly smaller then in inland vehicles 
or aircrafts, 

- the problem of the fast stopping of the ship movement or 
changing its direction due to the small water friction, 

- the possible large forces and moments from wind and 
waves, comparable with forces and moments from ship’s 
thrusters. 


But the most important problem is related to the fact that 
the ship is strongly nonlinear, multivariable object with time 
varying parameters. It’ s properties depend on load condition, 
speed value and direction, drift angle, trim, water depth, 
proximity of quays, other vessels etc. Very hard cross-coupling 


between all inputs and outputs also can be observed i.e. thruster 
action can change surge, sway or yaw simultaneously. 

Therefore any mathematical model of the ship dynamics 
only can approximate the real behaviour of the vessel. 
Consequently the testing processes of the various control 
systems dedicated for the ship’s motion steering cannot be 
performed only by simulation runs. 

Particular, it refers to the simpler linear models e.g. transfer 
functions or state space ones. More complex, nonlinear models 
(so called simulation models) are more accurate but researching 
results obtained via such tools also have limited accuracy. 

One can state that the simulation testing should be the first 
stage in the valuation process of any control systems but not 
the last one. Only wide range tests with the real ship during 
real-time experiments can guarantee that proposed steering 
algorithms are proper and likelihood. 

It was the main reason to decide in Ship Automation 
Department about building a floating model of the real 
vessel. 

Similar although bigger ships models are exploited in 
the Foundation for Safety of Navigation and Environment 
Protection at the Silm lake near Ilawa in Poland. Researching 
team from Ship Automation Department performed many tests 
with mentioned models in the past. Their results were published 
in e.g. [1, 2, 5]. But ship models exploited on the Silm lake 
are designated for improvement of captains and chief officers 
skills related to perform manual manoeuvres in different 
constrained areas but not for scientific researching. Therefore 
they are not equipped in a few navigational instruments e.g. 
radar, electronic maps or AIS. The researching ship described 
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here will be provided in needed devices as it will be shown in 
a further section. 

Constructed vessel is designated in the future to two main 
types of tests. The first one is related to different algorithms 
of the ship’s motion steering e.g. trajectory tracking, precise 
movement with slow speed and any drift angle, parallel motion 
of two vessel etc. Anti-collision algorithms, and agent systems 
state the second type of provided tests. 

A few persons from Ship Automation Department 
besides authors take part in some construction works with 
the researching ship at the beginning: Krzysztof Dziedzicki, 
Andrzej Januszewski, Roman Smierzchalski and Marcin 
Tobiasz. 


THE DESCRIPTION OF THE 
RESEARCHING SHIP - MANEUVERING 
SYSTEM 


The hull of this ship was built in Ship Design and Research 
Centre — Ship Hydrodynamics Division as a one from big 
number of hulls models of the merchant ships. 

It was constructed by means of thick horizontal slices of the 
wood, especially glued one slice to the another. The silhouette 
of the wooden hull is presented in Fig. 1. 


1.04|m 


5.80 m 


Fig. 1. The silhouette of the wooden hull 


Such hull was designated mainly to the water drag tests 
in Towing Tank Centre in Gdańsk-Oliwa. Therefore it was 
not suitable (in this pure wooden form) to installation of the 
thrusters and navigational devices. 

The transformation process of the hull, enabled the creation 
of the researching ship which can afloat on the lakes and the 
sea proceeded in a few steps. 

The first step was related to the cleaning of the hull, 
removing old pain and planning the outer side of the hull. 

In the next step the steel cage was constructed. The outward 
shape of the cage was accurately fitted to the inward shape of 
the wooden hull. All was joined together by screws (Fig. 2). 


Fig. 2. The steel cage inside the hull. Notice the construction 
of the hull by red painted slices of wood 


12 POLISH MARITIME RESEARCH, Special Issue 2012/S1 


The third step was related to the two holes made in the 
underwater part of the hull: first hole in the bow and the second 
one in the stern. They were assigned for tunnel thrusters. 
The next hole at the end of the hull was designated for main 
“engine” screw. 

Such prepared hull was covered by epoxide resin with a fibre 
glass mat and white painting special polyurethane paint. 

The general intention of the builders was likening of 
the constructed researching ship to the real ship named m/s 
“Finntrader” from Gdansk Shipyard. The side view id the 
mentioned ship is presented in Fig. 3. 

It should be emphasized that the researching ship is not 
a physical model of m/s “Finntrader”, although the ratios of the 
lengths and the beams are similar. Approximate scale is 1:30. 


Fig. 3. The silhouette of m/s ‘‘Finntrader” from Gdansk Shipyard 


The planned silhouette of the researching ship with the 
superstructure and hath covers is shown in Fig. 4. 


Fig. 4. The silhouette of the researching ship 
with the superstructure and hath covers 


The superstructure has the open construction (only the front 
part and the both side ones) due to necessity of the creation 
of the stands for two persons. There is also the possibility 
to mount, for particular tests a short lattice mast for radar 
antenna. The arrangement of the crew positions and the mast 
are presented in Fig. 5. 


Fig. 5. The side view of the researching ship with the crew stands. 
The removable mast are drawn by dashed line 


Described ship is equipped with battery-fed electrical 
thrusters: one main engine and two tunnel thrusters. 

The main engine for ship was adopted from standard 
outboard electrical drive for touristic boats made by Torqeedo 
(Fig. 6) [10 ]. The useless upper part of the drive was cut out. 

The specifications of the drives are as follows: 

a) electrical motor 


Base Travel 801L 


The researching ship “Gdynia“ 


-E 
A 


Fig. 6. The Torqeedo electrical drive Base Travel 800L [10] 
b) main propeller — diameter 30.48 mm 


The arrangement of the main drive elements are presented 
in Fig. 7 [4]. 
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Fig. 7. The elements of the main drive in the hull [4] 


The electric telegraph enables manual or computer steering 
of the main engine (Fig. 8). 


telegraph 
potentiometer 


manual 
switch 


besat 


PC 
Control System 
Fig. 8. Electric diagram of the main engine telegraph 


Tunnel thrusters: 
a) electrical motor 


Base Travel 801L 


b) thrusters screw propeller — type V-SET0086 


The block diagram of the tunnel thrusters electric system 
is presented in Fig. 9. 
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Fig. 9. The block diagram of the tunnel thrusters steering system 
The following elements was used in the tunnel thrusters 


electric system (see Fig. 9): 
a) mosfet type IRFP4468PbF from International Rectifier 


b) mosfet driver type TC4420 from Microchip 
c) encoder type MOK30 from P.P.H Wobit 


The microcontroller and other electronic elements were 
collected in one cover for every thrusters as it is presented in 
Fig. 10. 


Fig. 10. Electronic elements for tunnel thrusters drive 


The ship is also equipped with blade rudder made from 
stainless steel. The rudder is driven by Vetus hydraulic system 
3, 11] consists of: 

hydraulic pump type V-HTP2008R, 

- electrohydraulic pump typeV_EHPA12R2, 
- cylinder (servo-motor) type V-MTC3008, 
- nylon tubes, 

- rudder angle sensor type RFUI718. 
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The block diagram of the rudder driven system is presented 
in Fig. 11 and its installation arrangement photo in Fig. 12. 


Fig. 11. The Vetus blade rudder driving system installed on researching 
ship: 1) cylinder, 2) hydraulic pump, 3) nylon tubes, 4) electrohydraulic 
pump, 5) double check valve [11] 
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Fig. 12. The arrangement of the rudder driving system elements inside the 
researching ship. Notice the fragment of the steel cage [3] 


THE DESCRIPTION OF THE 
RESEARCHING SHIP — NAVIGATION 
SYSTEM 


To verify the ship’s control system one can propose to use 
areal model, equipped with real navigation devices. Model 


tests should reflect the real phenomena in a manner as close as 
possible to reality. In other words, modelling the phenomenon 
should be modelled similar to the phenomena taking place 
in reality. Therefore, building a physical model modelled 
on the real object, consider the need to preserve geometric 
similarity, kinematic and dynamic between the phenomenon 
being modelled and the modelling [7].These laws of similarity 
should be considered during the implementation of physical 
model tests of a moving ship in a real environment. Therefore, 
it is an important element to ensure proper operation of all 
equipment and components used to move the ship as well 
as to navigate. The manoeuvrability of the model are to 
provide: electric main drive, two bow thrusters and a hydraulic 
machine control. To observe the situation of the navigation on 
the sea area is conducted by: a system of radar, AIS, ARPA, 
electronic charts, GPS, log, echo sounder, anemometer, and 
VHF communications system. All the collected measurement 
parameters are transmitted by the card drivers and USB 
ports to a central control system through which it can be 
verified developed algorithm or control system. An additional 
advantage of using this model is the ability of the physical 
presence of researchers on board. This allows the continuous 
supervision and control of the working systems and also the 
fast identification of any irregularities. 

The structure of the measuring system used for the 
verification process, the ship’s control system is presented in 
Fig. 13. 

Developed system to verify the ship’s control system 
includes a number of devices associated with the conditions 
prevailing on the hydro meteorological sea area as well as 
equipment and actuators associated with control of the actual 
model. The model uses real navigation devices which are 
used on the bridges of merchant vessels. Mainly FURUNO 
equipment were used here (see Fig. 14) [8]. 

To determine the navigational situation prevailing at sea 
area measuring Furuno system was used. The system includes 
following devices: Marine Radar Unit RPU-015 with electronic 
maps and ARPA, AIS Transponder FA-150 and Hemisphere GPS. 
a variety of navigational information like own ship status, other 
ships status, radar plotting data, wind, water temperature, depth 
and other information from sensor are given by the system. 

In Marine Radar Unit RPU-015 the target detection is 
realized by sophisticated signal processing techniques such as 
multi-level quantization (MLQ), echo stretch, echo average, 
and radar interference rejecter. It is possible to adjust two guard 


Anemometer 


A 


Fig. 13. Configuration of the devices installed on board 
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Fig. 14. Configuration of the navigation system devices 


zones at required ranges in any sector. Other ship’s movements 
are assessed by advanced target tracking software and alerted 
by CPA/TCPA data readouts. Operation of the radar system 
at X-band and S-band allows the acquisition of echoes with 
an accuracy of 1% of range in use or 10 m whichever is the 
greater. 

The using of the maps allows to visualize the situation 
around the ship by plotting coastline, own ship safety contour, 
isolated underwater dangers, buoys, traffic routing systems, 
prohibited areas, fairways as required by IMO. A map is 
a combination of map lines and symbols to aid route planning 
and monitoring on the radar equipment. 30 nav lines may be 
stored and each line may contain up to 30 waypoints. Five 
nav lines may be simultaneously shown on the display. 200 
waypoints are available. Own ship and other ship tracks may 
be stored at a selected interval for repeated use. 

The moving objects are acquired automatically by means of 
ARPA subsystem. View the traffic parameters for each of the 
objects is possible by selecting it on the screen. The operator can 
select and display on-screen performance for the six objects at 
once. ARPA also allows: acquisition 100 targets automatically 
or manually; automatic tracking of all acquired targets on 
the display in the range from 0.1nm to 32 nm; marking from 
5 to 10 past positions of tracked targets at intervals of 30 s, 
1; 2; 3; 6 min; define suppression areas by combined with two 
acquisition areas of 3-3.5 and 5.5-6 nm, or 0.5 nm deep sector 
or circle in 0.3-32 nm; visualization vector true or relative 30 s, 
1; 3; 6; 15; 30 min for prediction of target motion; setting the 
collision warning CPA limit: 0.2-10 nm, TCPA limit: 0-99 min; 
define guard zone by sector or polygon who may be set in any 
effective area; and trial manoeuvre realized by dynamic or 
static, with selected delay time. 

The AIS device allows to obtain information about others 
ships - their type, charge, size and motion parameters. The 
system also allows the differentiation of individuals: sleeping, 
activated, dangerous, selected, lost targets. At the same time it 
is possible to collect information to 1,000 units in the range of 
VHF system to the sea area of certain frequencies to refresh 
the data in the range from 30 s to 60 min. AIS allows to define 
criteria of call CPA or TCPA alarm [8]. The Furuno AIS front 
panel is presented in Fig. 15. 


Fig. 15. The AIS front panel 


Measurement of the model speed and course is implemented 
by the GPS system and a log. The measurement system of 
Hemisphere company ensures a 2D position fix heading 
accuracy better than 0.1 degree rms with differential 


Fig. 16. The elements of the GPS device [9] 
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positioning accuracy of less than 60 cm, 95% of the time. 
Integrated gyro and tilt sensor deliver fast start-up times and 
heading provide updates during temporary loss of GPS. Fast 
heading and positioning output rates up to 20 Hz. Differential 
options including SBAS (WAAS, EGNOS, etc.) and beacon 
differentia. COAST ™ technology maintains accurate 
solutions for 40 minutes or more after loss of differential 
signal (Fig. 16) [9]. 


The WINDOBSERVER II anemometer from GILL 
company is used to determine the velocity and direction of 
the wind acting on the researching ship model (Fig. 17.), This 
device is characterized by the following parameters [12]: 

- measuring range for the wind: 0+ 65 m/s, 

- accuracy of measurement for wind speed: 2%, 

- measuring range for the wind direction: 0 + 359°, 

- accuracy of measurement for the wind direction: 2%, 

- operating temperature: -55sC to +70sC (with heated 
option). 


Fig. 17. The WINDOBSERVER II anemometer antenna [12] 
The PC Control System 


The modern software installed on a PC included the 
measuring subsystem and also the control one creates the full 
manoeuvring system. The test, measuring and control signals 
via the USB and card transmitters are supplied for all devices. 
The signals are transmitted in the standard NMEA 0183 and 
standard + 5V. In addition to transferring data to and from the 
model, you can use GPRS. Other information signals such as 
signals from the log, echo sounder and anemometer are fed via 
RS422 or link directly to the navigation system Furuno. The 
computer control system allows also to connect signals from 
other devices, for example wave meter or thermometer, etc., 
to determine the conditions aroun the ship. 

All information can be presented on the screen inside the 
ship (Fig. 18). 


PLANNED TESTS 


As it was mentioned in the first section of the paper the 
researching ship is designated to two main types of tests: 
1. The different algorithms for the ship’s motion steering will 
be checked taking into account the influence of the wind 
and waves. 
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Fig. 18. View of the exemplary navigation situation 
from PC Control System [6] 


2. Anti-collision algorithms, and agent systems state the 
second type of provided tests and will used presented 
navigational instrumentation (radar, ARPA, electronic 
maps, AIS, etc.). 


The computer algorithms enables full ship motion control 
can be associated with following control subsystems: 
a) heading and/or speed stabilization, 
b) stabilization of the turning radius, 
c) roll damping, 
d) trajectory tracking, 
e) precise steering with the low speed, 
f) dynamic ship positioning (DSP), 
g) automatic single-buoy mooring, 
h) turret-mooring for FSO, FPSO or FPDSO ships. 


Systems denoted by letters a) + d) can be treated as 
one-dimensional control units, but the others only as 
multidimensional ones. The heading, speed, turning radius 
stabilization and trajectory tracking can be performed by means 
of blade rudder or main propeller. Precise steering and DSP 
processes needs the ship equipped with at least a few driving 
devices on bow and stern. 

The set of the driving devices installed on researching ship 
enables using any ship motion control subsystems. 

There are planned among others the following tests: 

a) the heading stabilization by means of adaptive, neural and 
fuzzy logic autopilots, 

b) the trajectory stabilization by means of indirect “line of 
sight” controller, 

c) the translation of the vessel with low speed and any drift 
angle by means of multidimensional regulator based on 
robust control theory or LMI approach. 


Planned exemplary manoeuvres are presented in Figs 19 
+ 21 below. 

The second type of research is dedicated for verification 
and develop a reliable anti-collision algorithms and expert 
systems for navigators. Modern maritime navigation requires 
the navigator to determine the optimal speed and routes, also in 
danger of conflict. It is not always the task is completed within 
a reasonable degree as evidenced by the statistics of accidents 
at sea. Simultaneously the problem of determining optimum 
routes for vessel between its current position and given final 
destination point, including the situation occurring when 
navigating a particular route to overcome, was the subject of 
many studies and has not been finally resolved. To determine 
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the optimal (economical and safe) routes and steering the vessel, 
it is proposed to test the many control algorithms including 
agent systems. The use of the agent system of navigation is 
a new approach to the implementation of collision avoidance 


systems at sea. 
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Fig. 20. The changing of the quay with rotating manoeuvre 
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Fig. 21. The leaving harbour manoeuvre 


The agent system is the one that consists of a number 
of agents that cooperate with each other in order to solve 
a particular problem in this case - maximum level of safety 


at sea. The proposed system would automatically pursue 
three fundamental tasks related to ensuring safety at sea: data 
collection and analysis of the current situation around the ship’s 
own navigation, automatic negotiation between ships in the sea 
area or shore stations to determine the potential area in which 
it would be the route transition (the job has not yet found in 
decision support systems) define routes and its auto correction 
depending on changes in current navigational situation. The 
Particular platforms of agent system Could be installed at 
different ships or waterside stations [6]. The agents work 
within a single platform of the agent system for exchanging 
data between them (Fig. 22). 


Fig. 22. The structure of the agent system 


Practical implementation could have promising effects in 
the form of relief work of the navigator, increase safety at sea, 
and reduce operating costs of the vessel. Agent system in its 
simplicity, offers the opportunity to review other methods and 
algorithms to determine optimal routes including anti-collision 
manoeuvre. 


CONCLUSIONS 


The proposed process for reviewing the ship’s control 
system is an indirect method of verification between the 
low-cost methods of conducting research using computer 
simulations and tests with the use of expensive full-sized ships. 
Despite the use of the model constructed on a large scale, 
potentially results obtained may be considerably closer to the 
phenomena occurring in reality than the results that may be 
obtained through computer simulations. Thanks to electronic 
navigation modelling environment, which moves in the actual 
model of the ship, it is possible to very accurately modelling 
a wide variety of situations that contain various combinations 
of navigational limitations of a static (land, shoals, channels, 
fairways, engineering structures, etc.) and dynamic (icebergs, 
the other moving ships and facilities, etc.). Application 
submitted verification process decision support systems and 
vessel traffic management can contribute significantly to the 
improvement of maritime safety and reduce operating costs 
of the ships. 
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Decision support system 
for collision avoidance at sea 


Agnieszka Lazarowska, M.Sc 
Gdynia Maritime University 


ABSTRACT 


The paper presents design and realization of computer decision support system in collision 
situations of passage with greater quantity of met objects. The system was implemented 
into the real ship electro-navigational system onboard research and training ship m/v 
HORYZONT II. The radar system with Automatic Radar Plotting Aid constitutes a source 
of input data for algorithm determining safe trajectory of a ship. The article introduces 
radar data transmission details. The dynamic programming algorithm is used for the 
determination of safe optimal trajectory of own ship. The system enables navigational 
data transmission from radar system and automatic determining of safe manoeuvre or safe 


trajectory of a ship. Further development of navigator 5 decision support system is also presented. Path 
Planning Subsystem is proposed for the determination of global optimal route between harbours with the 
use of Ant Colony Optimization algorithms. 


Keywords: anti-collision system; collision avoidance; safety of navigation; safe ship control; computer simulation; 
decision support system; navigational data; marine navigation; path planning; ant colony optimization 


INTRODUCTION 


Technological development led to an increased marine 
traffic, which caused navigation to become more demanding for 
deck officers. Assertion of safe navigation simultaneously with 
minimal operational costs constitutes the main issue in present- 
day marine transport. Carrying out navigation represents 
a complex process, because it requires continuous analyzing 
of huge amount of information. Incorrect assessment of current 
navigational situation can lead to collision situations often with 
very tragic consequences. Therefore it is necessary to support 
deck officers in anti-collision decision making process. 

Currently onboard a ship the function of anti-collision 
system is accomplished by Automatic Radar Plotting Aid 
(ARPA). According to the requirements of Safety of Life at 
Sea Convention (SOLAS), enacted by International Maritime 
Organization (IMO), all ships of 10 000 gross tonnage and over 
constructed on and after 1 July 2002 have to be equipped with 
Automatic Radar Plotting Aid [1]. ARPA system provides the 
possibility to track automatically at least 20 targets. The system 
also generates dangerous target alarm which indicates that the 
computed values of Time to the Closest Point of Approach 
(TCPA) and Distance of the Closest Point of Approach (DCPA) 
exceed the specified safe limits. ARPA system offers also 
a collision avoidance support function called trial manoeuvre. 
By using this function the deck officer has the opportunity to 
check the effects of own ship planned manoeuvre. It enables 
simulation of course change as well as speed change of own 


vessel [2]. Therefore there is still a possibility to improve 
performance of anti-collision system by implementation of 
a decision support system automatically determining the safe 
trajectory of own ship. 

Human error causes 75-96 % of marine accidents [3]. 
Application of a system eliminating human subjectivism 
in decision making process will cause the number of ship 
accidents to decrease, providing protection of human life and 
health, transported cargo and natural environment. Moreover 
the system will also assure economical transport by taking into 
account optimality criterion in the form of the smallest lost of 
way on overtaking other ships. 


DECISION SUPPORT SYSTEM 


The role of the decision support system is to aid the deck 
officer in the process of choosing the anti-collision manoeuvre. 
The main function of the system is to automatically determine 
safe manoeuvre or safe trajectory of own ship. It is achieved 
by application of computer algorithms. In the designed system 
input data for the algorithm are collected from Automatic Radar 
Plotting Aid. 

Designed Navigator’s Decision Support System (NDSS) 
is composed of two subsystems - Navigational Situation 
Display System (NSDS) and Determining Safe Trajectory 
System (DSTS) (Fig. 1). The first one collects input data from 
ARPA system and displays current navigational situation. 
It also contains a system log to save the data received from 


POLISH MARITIME RESEARCH, Special Issue 2012/S1 19 


Decision support system for collision avoidance at sea 


ARPA system. The second subsystem is responsible for the 
performance of anti-collision calculations. Initialisation 
of this subsystem is achieved manually by the operator or 
automatically, after activation of dangerous target alarm in 
ARPA system. 

Foundations of the system include collecting of input data 
describing current navigational situation from Automatic Radar 
Plotting Aid, taking into account specific character of safe 
ship control process, characterized by great course changes in 
range from 20° to 90° and reduction of speed of no more than 
30% and meeting requirements of International Regulations 
for Preventing Collisions at Sea (COLREGS). It also takes 
into consideration weather conditions distinguishing between 
situations of good and poor visibility and includes dynamic 
characteristics of own ship [4, 5]. 

To sum up, the system provides automatic anti-collision 
decision support in near real-time, reliable operation, not 
causing disruption in operation of basic anti-collision 
system, user-friendly and clear presentation of proposed 
manoeuvre. 

Tasks realized by the system comply automatic collection 
of data concerning current navigational situation, analysis of 
current navigational situation, indicating dangerous situations 
and showing present level of navigational safety and automatic 
determination of anti-collision manoeuvre. 


p, [RADAR | Loc | GYROCOMPASS| 


a TE | 
ARPA (Automatic Radar Plotting Aid) | 


v, ¥, Vj, ¥j, Nj, Dj 


NDSS — Navigator’s Decision Support System 


NSDS — Navigational Situation Display System 


DANGEROUS TARGET en 


aK fax. AV, Ps(A¥e, AVe, t) 
operator $ action 


DSTS — Determining Safe Trajectory System 


Fig. 1. Block diagram of decision support system. Ds — safe distance, 
Y — own ship course, V — own ship speed, Pj — tracked target course, 
Vj — tracked target speed, Nj — tracked target bearing from own ship, 
Dj — tracked target distance from own ship, A¥ — manoeuvre of course 
change, AV — manoeuvre of speed change, Ps(A Ve, AVe, t) — safe ship 
trajectory 


AUTOMATIC RADAR PLOTTING AID 
AS A SOURCE OF DATA FOR ANTI- 
COLLISION CALCULATIONS 


On the basis of signals transmitted from radar’s transceiver 
extractor determines position of detected targets - bearing N; 
and distance D, from own ship. In the next step, computer 
system calculates motion parameters - course ¥. and speed V. 
of tracked targets and approach parameters - distance of the 
closest point of approach DCPA and time to the closest point of 
approach TCPA. Despite the signals from radar system ARPA 
receives also data from gyrocompass and log about own ship 
course and speed. 

To summarize, ARPA system constitutes a source of the 
following parameters: course and speed of own ship, target 
vessel course, speed, bearing and distance from own ship 
(Fig. 2). These parameters compose input data to navigator’s 
decision support system. Despite these parameters, in the input 
to decision support system there has to be included information 
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about safe distance Ds. This parameter is estimated by the 
navigator and describes minimal acceptable distance of passing 
met objects by own ship. 


RADAR 


LOG ||GYROCOMPASS 
J 


TRANSCEIVER 


CONTROL PANEL 


POWER SUPPLY 
UNIT 


Fig. 2. Block diagram of radar and ARPA. ¥ — own ship course, 

V — own ship speed, Yj — tracked target course, Vj — tracked target speed, 
Nj — tracked target bearing from own ship, Dj — tracked target distance from 
own ship, DCPA — tracked target Distance of the Closest Point of Approach 

TCPA - tracked target Time to the Closest Point of Approach 


SPECIFICATION OF RADAR DATA 
TRANSMISSION TO DECISION SUPPORT 
SYSTEM 


The data transmission is performed in accordance with IEC- 
61162-1 standard specified by International Electrotechnical 
Commission and is in a close convergence with NMEA 0183 
standard. The IEC-61162-1 standard defines an interface 
requirements for data communication between maritime 
electronic instrument and external device. Communication is 
performed with the use of serial asynchronous data transmission 
of RS 232 type. Applied baud rate is 4800 bits per second. The 
data frame consists of one start bit, 8 data bits and one stop 
bit. There is no parity bit. The data frame structure is shown 
in Figure 3. Besides the definition of data frame format, IEC 
61162-1 standard characterizes also the sentence structure. 
Figure 4 presents an example of NMEA data received from 
GPS eTrex Legend. 

It is possible to transmit different navigational data, among 
others, concerning geographic position of our ship - GLL, 
depth below transducer - DBT, wind direction and speed 
- MWD or rate of turn - ROT. The selection of transmitted 
sentences is configurable and can be changed with the use of 
service menu of our ARPA system. Via the service menu we 
can also configure the transmission frequency of every output 
sentence. In SAM Electronics Radarpilot 1100 system it can be 
selected from 4 different values: 1,5, 10 or 30 seconds, when 
the standard repetition time of OSD sentence is 1 second and 
of TTM sentence is 10 seconds or more. 

Data needed for anti-collision calculations are included in 
sentences with sentence formatters marked OSD and TTM. 
The first one describes own ship data when the second one 
defines tracked target message. In the sentence before the 
sentence formatter appears a talker identifier that specifies 
the device which constitutes the source of data. When the data 
are received from radar system including ARPA, the talker 
identifier is marked with RA. In sentences defining own ship 
and tracked target data there can also be found characters used 
for marking reference system on which the calculations of 
vessel course and speed are based, B means bottom tracking 
log, M stands for manually entered, W for water referenced, 
R means radar tracking and P positioning system ground 
reference. There are also symbols representing speed units, 
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K means km/h, N stands for knots and S for statute miles/h. 
In addition tracked object message defines target status, which 
can be marked as L for target that has been lost, Q for query, 
which means that the target is in the process of acquisition 
and T for tracking. Moreover, tracked target message includes 
information presenting type of acquisition. It can be marked 
with a for automatic, M for manual or R for reported, for 
example from AIS [6, 7, 8]. 


MEI MEI Rbk 


Start Stop 
bit Data bits bit 


Fig. 3. Data frame structure 


The sentence structures concerning own ship data and 
tracked target message are as follows: 


$RAOSD,a.a,A,c.c,X,X.X, Y,Y.Y,Z.Z,N*hh<CR><LF> 


RA - talker identifier: RA - radar 

OSD - sentence formatter: OSD - own ship data 
a.a - heading, degrees true 

A - heading status: a -valid, V -invalid 

c.c - vessel course, degrees true 

X - course reference: W -water-referenced 
x.x - vessel speed 

Y - speed reference: W -water-referenced 
y.y - vessel set, degrees true 

z.z - vessel drift (speed) 

N- speed units: N - knots 


$RATTM,aa,b.b,c.c,T,t.t,x.x, T,y.y,z.z,X,n, Y,R,hhmmss. 
ss,M*hh<CR><LF> 


RA - talker identifier: RA - radar 

TTM - sentence formatter: TTM -tracked target message 

aa - target number: 00 to 99 

b.b - target distance from own ship 

c.c - bearing from own ship 

T - degrees true 

t.t - target speed 

x.x - target course, 

y.y - DCPA - distance of closest point of approach 

z.z - TCPA - time to closest point of approach [min] 

X - speed, distance units 

n- user data (eg. target name) 

Y - target status: L-lost, Q - query (target in the process of 
acquisition), T -tracking 

R - reference target (used to determine own ship speed) 
hhmmss.ss- time of data (UTC) 

M - type of acquisition: a -automatic, M-manual, R-reported 
(eg.AIS) 


-6,250112,3.7, E A*1A 
18 478 ,.259.4, .U. AGA 
3.8,4,34.6,M,,*?7D 


E 
-86,17 00. ig” "24, 298.00.19.16, 046, 4076 
Pigesae. A, Ax4P 


CPLL? 541 BÉ8A,N, 18 
GPBOD,207.3,T,203.6,M, SAM, 
GPUTG. 295 .6.1.291.9,M.0.8.N,.8.8,K™45 


Fig. 4. NMEA data received from GPS eTrex Legend produced by Garmin 
with the use of application implemented in C programming language 


COMPUTER PROGRAMME FOR 
COLLECTING RADAR DATA AND 
CALCULATING SAFE TRAJECTORY 
OF OWN SHIP 


Various methods were proposed by many researchers for 
the determination of anti-collision manoeuvres, among others, 
game control methods [9, 10, 11, 12, 13], neural networks [14], 
fuzzy logic [15, 16, 17, 18], evolutionary and genetic algorithms 
[19, 20] were applied. 

Primary implementation of the designed system was 
conducted with the use of dynamic programming algorithms 
for determination of safe optimal trajectory of own ship 
implemented in Matlab environment [21, 22]. Applied dynamic 
programming method takes into consideration dynamics of 
own ship and requirements of International Regulations for 
Preventing Collisions at Sea, what is essential in safe ship 
control process. The application executing transmission of 
anti-collision data was also designed in Matlab software and 
integrated with the algorithm determining the safe trajectory 
of own ship. 

First two steps of communication algorithm are formed by 
the definition of serial port used for communication with ARPA 
system and configuration of its parameters. Transmission speed 
described as baud rate has to be set to 4800 bits per second. 
The next action to be taken is opening the connection. From 
that moment all available data are received via serial interface. 
If it is not possible to restrict transmitted sentences only to 
these with OSD and TTM sentence formatter, because access 
to service menu of radar system is unavailable, there exists 
a need to filter the data stream to obtain only strings marked 
with ‘$RAOSD’ and ‘$RATTM”’. These obtained strings have 
to be parsed then to select data used as an input to anti-collision 
calculations. So parsing routines have to be used to obtain 
own ship course and speed from OSD string and target vessel 
course, speed , bearing and distance from own ship from TTM 
string. After acquisition of necessary data the connection can be 
closed and action passed to anti-collision algorithm, as a result 
of which consecutive course and/or speed changes of our ship 
are obtained. Consecutive steps executed by communication 
algorithm are shown in Figure 5 [23,24]. 


EXPANDED DECISION SUPPORT SYSTEM 


Harris, Tran and Wilson proposed a collision avoidance 
system called Maritime Avoidance Navigation, Totally 
Integrated System (MANTIS), which is composed of the 
Localisation of Vessel States and it’s Environment (LVSE), 
Automatic Collision Avoidance Advisory Service (ACAAS), an 
Integrated Display System (IDS), Path Planning and Scheduling 
Service (PPSS) and Automated Ship Guidance and Control 
(ASGC). Their proposal stimulated further development of 
decision support system presented in this paper [25]. 

Decision support in the ship control process should be 
composed of two tasks. Besides advising in anti-collision 
decision making process, it should also constitute a tool for 
determination of optimal route of passage between harbours. 
For that reason an extended version of decision support 
system is proposed. The system of hierarchical, multilayer 
structure includes Path Planning Subsystem (PPS) for the 
determination of global route from departure to destination 
harbour and Determining Safe Trajectory Subsystem (DSTS) 
for the determination of local route. DSTS system operation 
is initialized manually by the operator or automatically, 
when dangerous target alarm is activated in ARPA system. 
It calculates anti-collision trajectory from current position 
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setting simulation parameters 


creating a serial port object 


setting transmission parameters 


opening connection 


collecting own ship data 


parsing the string to obtain ¥ and V 


z 
z 


collecting tracked objects data 


parsing the string to obtain ‘Yj, Vj, Nj and Dj 


closing port 


anti-collision calculations 


displaying results 


Fig. 5. Diagram of communication algorithm integrated with decision 
support system. Y — own ship course, V — own ship speed, Yj — tracked 
object course, Vj — tracked object speed, Nj — tracked object bearing from 
own ship, Dj — tracked object distance from own ship 


of own ship to the next waypoint determined by global path 
planner. This local route is determined with the use of dynamic 
programming algorithm. While for the determination of global 
path in PSS algorithm Ant Colony Optimization is proposed. 
Global path planning algorithm takes into account navigational 
obstacles in the form of land, canals, shallows etc. 

Ant Colony Optimization operation principle is based on 
ant colonies foraging behaviour. It has been developed that 
ants, while moving between food source and their nest, deposit 
chemicals on the ground called pheromones. Other ants smell 
this pheromone trail and choose the path marked by stronger 
pheromone concentrations. Investigation of this trail-laying and 
trail-following behaviour of some ant species have shown that 
this indirect communication called stigmergy forms a method 
in which ants find the shortest path between food source and 
their nest. 

ACO algorithms belong to the branch of artificial intelligence 
called swarm intelligence. Swarm intelligence is concerned 
about design of multi-agent systems. It takes inspiration from 
the collective behaviour of social insects such as ants, bees, 
wasps and termites. In analogy to real ants behaviour, ACO is 
based on indirect communication of a colony of simple agents 
called artificial ants. The ACO approach to solving optimization 
problems was first introduced by Dorigo in collaboration with 
Colorni and Maniezzo in the early 1990’s. 
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The representation of a problem with the use of ACO is as 
follows. If a minimization problem (S, f, Q) is considered, S 
is the set of candidate solutions, fis the objective function and 
Q is the set of constraints. The objective function f assigns an 
objective function cost value f(s, t) to each candidate solution 
s € S. The aim is to find a globally optimal solution S, € S, 
it means a minimum cost solution that satisfies the constraints 
Q. The description of the problem exploited by the ants include 
a finite set C = {c}, Cz, ..., Cy Sof components. The problem 
states are defined as a sequence x = (c Cig very Cys ...) over the 
elements of C. X is a set of all possible sequences and |x| is 
the number of component in the sequence x, denoted as the 
length of a sequence x. The set of feasible states X with c X c 
X is defined by the finite set of constraints Q. A set of feasible 
solutions is denoted by S* (with S* c X and S* c S). Each 
candidate solution s € S have a cost f(s, t) associated. 

It spite of above-mentioned definition of the problem, it 
can be considered, whether to implement constraints Q in 
a hard way allowing ants to build only feasible solutions or 
in a soft way, where ants can also build infeasible solutions, 
which will be penalized then, in dependence of their degree 
of infeasibility. 


The travelling salesman problem was the first application of 
ACO algorithm. In this path optimization problem the goal is to 
find a minimal length closed route connecting m given cities, 
on the assumption that each city can be visited only once. 

The problem definition can be presented in Euclidean space, 
where the distance between two cities i and j is described as 


follows: 
dis f(x; zx +: =y (1) 
where: 


x; and y;— the coordinates of city i, 
x; and yj— the coordinates of city j. 


It can also be described as a graph (N,E), where the nodes 
N are defined by the cities and the connections between the 
cities are presented as the edged E of the graph. 

Given the problem representation as stated above, during 
every iteration of algorithm, each ant k (k = 1, ..., m) moves 
on the problem graph from one city to another by applying 
a probabilistic transition rule, until it completes the tour. The 
transition rule is expressed by: 


lee) 
miO- nf 
o0 
Den EO nt 
where T(t) is the amount of virtual pheromone trail on the 
edge connecting city i to city j, which reflects the experience 
acquired by ants during problem solving, nj= I/d; is a heuristic 
information called visibility. Parameters a and B determine the 
relative influence of pheromone trail and heuristic information 
and Nk; is the set of cities which ant k has not yet visited, when 
it is in city i. When the solution construction at present iteration 


ends, the pheromone trails are updated. The expression defining 
this procedure is given by: 


tut D=C-p) ty) + D At G) 


pi (t)= ifjeN, © 


where p is the pheromone trail evaporation rate, which enables 
the algorithm to “forget” bad decisions previously made by the 
ants, At k (t)is the amount of pheromone ant k deposits on arcs 
belonging to its tour and is defined as: 
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ae if arc(i, j)is used by ant k 


At) (t) = (4) 
0 othervise 


Ship’s path planning problem is represented by a graph (N, 
E), where nodes N are placed in acceptable waypoints and edges 
E connecting these nodes, for which passing of a ship is possible 
without transgressing constraints [26, 27, 28, 29, 30]. 


CONCLUSIONS 


Designed decision support system enables to automatically 
determine safe manoeuvre or safe trajectory with the use of 
information from automatic radar plotting aid. The system 
provides possibility to implement different methods of 
determining safe ship trajectory in collision situation other than 
applied dynamic programming and to verify their performance 
in real navigational situation and in real electro-navigational 
system onboard a ship. Further research is conducted to 
make a system independent from Matlab environment and 
implemented in C programming language to constitute 
a stand-alone application (Fig. 8). Figs 6 and 7 present graphic 
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representation of navigational situation solved by C language 
implementation of dynamic programming algorithm. In an 
expanded version the decision support system could present 
safe trajectory on ARPA display or Electronic Chart Display 
and Information System (ECDIS) display with the proposed 
value of course and/or speed change of own vessel. Further 
development of the system could lead to full automation 
of the process of safe ship control. It could be achieved by 
sending rudder order after confirmation by the deck officer. 
Joined with the path planning system onboard the ship would 
enable automatic passage from point A to point B, because the 
system would be capable of automatically solving collision 
situations at sea. 
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Fig. 8. Results of ship s safe trajectory determination with the use of 
dynamic programming algorithm implemented in C programming language 
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ABSTRACT 


The paper presents software environement for fast 
prototyping and verification of motion control systems for 
ship. The environement is prepared for isomorphic reduced 
ship model which is used for training and in research in 
a area of ship motion control. The control system is build 
using Matlab-Simulink-xPC package which simplifies 
and accellerates design and verification of new control 


algorithms. The systems was prepared also for Hardware- 
in-the-loop trials when a designed control system is tested inside a virtual environment instead of real 
actuators, disturbances, communication and measurement devices. 


Keywords: ship controll Matlab; xPC Target; fast prototyping; hardware-in-the-loop; 
real-time control; computer simulation 


TRIALS ON A MATERIAL SHIP MODEL 


It is common rule that final verification of a designed control 
system should be performed in a real world on the real control 
object. Because of economical and safety reasons such trials 
of control systems for seagoing ships are difficult to perform. 
Unfortunatelly nowadays control systems engineers have to 
incorporate innovative technologies in short time which is 
limited by requirements of maritime industry. Therefore most 
of trials are made in computer simultion, when not all design 
factors can be tested and solved, so designed systems can not 
be fully robust and reliable. 

Interesting option for verification of designed control 
systems are tests with reduced material model of ship. Using 
rules of mechanical similarities results obtained on the reduced 
material model can be scaled to those which can be obtained 
on an original ship. Ship Handling, Research and Training 
Center on the Silm Lake in Iława, Poland has many reduced 
ship models in scale 1:24 and 1:16 which are mainly used for 
a captains training. But since they have been equipped with 
additional control devices these models are perfect control 
objects for research purposes. 

One of them called Blue Lady (Fig. 1) is used by research 
team from Department of Ship Automation at Gdynia Maritime 
Academy for verification of new control algorithms. Blue 
Lady is a reduced material model of the 176 000 DWT VLCC 
(Very Large Crude Carrier) tanker in a scale 1:24 made of an 
epoxy laminate [1, 2, 3]. It has propulsion system powered by 
electric motors and fed from set of batteries. The model was 


made subject to the rights of geometric, kinematic and dynamic 
similarity, so the results and measurement data obtained on the 
lake can be scaled to the results recorded on a real tanker. 


main tunnel thrusters 
propeller 


rotating pump thrusters 


Fig. 1. VLCC “Blue Lady” tanker model as a training ship 
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Tunel thruster H2 


Tunel thruster H3 


Stern rotating 
rudder thruster H4 


| cps NMEA 0183 / RS232 


NMEA 0183 / RS422 


NMEA 0183 / RS232 


49}|O4JUOD JAYS OL 
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Fig. 2. Control equipment of the training ship Blue Lady 


Fig. 2. presents devices used for automation control of 
the Blue Lady training ship: main engine with proppeler, 
conventional blade rudder and four variable speed jet thrusters: 
two of them are tunnel thrusters and other two are rotating 
thrusters with constrained turn angle. Direct control of actuators 
and measurement of real settings is performed by dedicated 
microcontroller system (called a hub) through set of 12-bit 
analog-digital and 10-bit digital-analog converters. The hub 
can cooperate with the master controller through bidirection 
interface RS-232. Data exchange format has form of binary 
records. 

Moreover a gyrocompass for heading measurement, 
a ultrasonic anemometer for wind disturbances, and DGPS 
(Differential GPS) receiver with a reference station fixed on 
a lake shore - for recording of actual position are installed on 
the training ship. These devices transmittes data in NMEAO0183 
format through RS232/422 interfaces 


CONTROL SOFTWARE 


For scientist and engineers in a field of automation control 
Mathwork’s MATLAB is well known and frequently used 
programming package. Matlab users can make calculations with 
its programming language with many toolboxes and also can 
simulate control systems behaviour with the powerful Simulink. 
Less known is knowledge that some toolboxes gives possibility 
to control systems in a real time. The toolboxes that are used 
for building of real control systems are Simulink Coder (former 
Real Time Workshop) [4] and xPC Target [5]. 

Simulink Coder creates ANSI/ISO C code from given 
Simulink diagram and Matlab functions. The code can be used 
for developing and testing purposes, particulary for acceleration 
of simulation, rapid protoyping, real-time applications and 
hardware-in-the-loop testing. User can tune and monitor the 
generated code using Simulink or run and interact with the 
code outside MATLAB and Simulink. The code created by 
Simulink Coder can be compiled to the executable for many 
target platforms. 
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Using xPC Target toolbox the code generated by Simulink 
Code can be executed on remote computer, which become 
a real time controller. xPC Target provides a real-time kernel 
and a library of drivers to include in a Simulink diagram. 
xPC Target supports a comprehensive selection of I/O which 
includes: analog and digital inputs/outputs in many formats, 
serial communication ports (RS232/422/485), incremental 
encoder, Ethernet (UDP/IP, real-time raw Ethernet), video 
cards (USB WebCam, CameraLink), communication data buses 
and protocols (e.g. CAN) and more. The toolbox additional 
features are a host-target interface for real-time monitoring, 
parameter tuning, and data logging. xPC Target computer 
can work under constant monitoring of an master computer 
running Matlab software but also it is possible to generate 
a stand-alone software which boots PC and then works without 
any supervision. 

Popularity of Matlab makes it optimal choose for preparation 
ofa control setup for a research purposes. Transition from pure 
computer simulation to real application is much easier, because 
the same Simulink diagrams used in simulation are used in 
generation of control software for final controllers. Moreover 
since performance of nowadays microprocessors still fast 
increase, even developing of commercial, industry control 
systems with Simulink-xPC Target toolbox becomes interesting 
option, because time of obtaining of a final solution becomes 
the most important factor in a design process. 

For building of xPC Target control systems two PC-type 
computers are required. One of them plays role of supervising 
computer (host) and runs MS Windows operation system 
and Matlab-Simulink software. Because of reqirements of 
nowadays Matlab version this PC have to be rather fast unit. 
Second PC computer (target) is used as a hardware platform 
which runs dedicated, propertiary real-time kernel and control 
software. Creator of Matlab - Mathworks presents xPC 
Target Turnkey programme which combines xPC Target with 
a high-performance real-time target machine and I/O modules. 
Target machines and I/O modules meet a wide range of I/O 
connectivity, form factor, and environmental requirements. 
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Also using of a standard low-cost computer is possible, but it 
is recomended to use a dedicated and suported by xPC Target 
toolbox extension cards, because writing of software drivers 
for non-standard extension cards requires much programmer 
experience. Target computer performance should be selected 
according to control algorithms requirements, but in many 
cases is not critical factor. 

Both computers: host and target have to be equipped with an 
Ethernet networks adapters and connected together via TCP/IP 
protocol. The network junction is used for control soffware 
transfer from host to target computer, and also for monitoring 
and supervising. It is possible to generate stand-alone software 
for target. In this configuration target computer boots real-time 
kernel and loads control software directly from a hard disk. 

Fig. 3 presents how Matlab can be used during design 
process of anew control algorithms. In a first stage of 
a design Simulink is powerful tool for off line data processing, 
identification of a mathematical model of a process and fast 
verification of designed controllers. In a last third stage the 
final verification is done during trials on a lake, where a PC 
controller is powered by a PC-Target real-time kernel which 


Matematical 
model 
of Blue Lady 
dynamics + 
interfaces 


executes controll algorithm defined by a Simulink diagram. As 
a graphical Simulink diagram can be simply modified therefore 
such software can be called a real fast-prototyping tool. 

Another speed-up is introduced by a middle second 
verification stage so called Hardware-in-the-Loop (HiL). 
Real control system consists not only of controler but also it 
consists of measurement systems, actuators, communication 
interfaces, etc. HiL is a methodology of testing a whole 
controller’s functionality and communication. HiL significantly 
reduce development time of control systems by testing the 
controller inside a virtual environment instead of real actuators, 
disturbances, communication and measurement devices, so all 
hardware and software components of control system can be 
tested before final trials. 


RESEARCH SETUP 


As shown on figure 2 training ship Blue Lady requires four 
RS232/422 serial ports for communication with supervision 
controller. Therefore a four-port multi I/O PCI card has been 
installed in the xPC Target computer. The card was selected 


Stage | 
Computer 
simulation 


Stage Il 
Hardware-in- 
the-loop trials 


PC Target 
Standalone mode 


Training ship Blue Lady 


PC Target 
Real-Time Kernel 


Stage Ill 
Verification 


4 x RS232/422 


ona 
reduced 
material 
model 


PC Target 
Real-Time Kernel 


Fig. 3. Verification of control algorithms for ships 
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from those recomended by Matlab — it is QSC-100-d9 four-port 

RS-232 Multiport Serial Adapter from Quatech. Converting 

between RS232 and RS422 interfaces is done by additional 

external converter. 

The xPC serial port driver works with streams of bytes 
which represent data received from serial port or transmitted to 
a port. Data streams from receiver must be decoded according 
to a proper data format. Gyrocompass, GPS and anemometer 
send data in a NMEA 0183 ASCII format which can be 
decoded using standard ‘ASCII Decode’ block, where format 
string is defined in a similar way like in function sprintf in a C 
language. 

Data received from the hub is packed into binary record 
which consists of: 

- two-byte header ‘BL’; 

- 14word-size numbers from A/D converters which measure 
actual control parameters, valid are only 12 less significant 
bits in every number; 

- 5 bytes from digital inputs; 

- control byte which is sum modulo 256 of all previous bytes 
in the record. 


The multi I/O card sends only only data to the hub in form 

of binary record which consists of: 

- three-byte header ‘MP’; 

- 8 word-size numbers for D/A converters which controls 
actuators (valid are only 10 less significant bits in every 
number); 


E] Sceleton4 


Fie Edt View Simulation Format Tools Help 
0 ees 


- 5 bytes for digital outputs; 
- control byte which is sum modulo 256 of all previous bytes 
in the record. 


For the purpose of trials on ship models universal template 
diagram in Simulink have been developed. The diagram is 
called a sceleton and is shown in Fig. 4. 

The template diagram contains the drivers for Multi I/O with 
four RS232 ports and boxes which represents data decoders and 
encoder (Figs 5, 6, 7). All data from input decoders are brought 
to the box Control Algorithm, where Simulink diagram which 
represents a verified algorithm will be included. Each team 
member has the ability to insert and modify his/her algorithm 
which is currently working on: to control the course of the 
vessel, its speed, movement on the trajectory and dynamic 
positioning, etc. and no others part of the sceleton have to 
be modified. A simple PID course controller is shown as an 
example of control algorithm on figure 8. Output from the 
Control Algorithm is brought to Output Encoder, which pack 
data that are send to the hub and controls ship actuators. 

The skeleton also provides monitoring of all decoder 
outputs, encoder inputs from target computer into a Matlab 
workspace in a host computer. This data can be processed by 
additional Matlab scripts and stored on your hard disk for later 
processing. 

Target scopes are components which enables local 
monitoring facilities for selected parameters on a screen of 
target xPC computer. 
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Fig. 4. Sceleton - Simulink control framework used in trials 
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HARD WARE-IN-THE-LOOP SIMULATION 


For Hardware-in-the-Loop (HiL) trials a virtual 
environement of system on the Silm lake has also been 
build using Matlab-Simulink-xPC Target tools. Devices and 
dynamics of areal traing ship are simulated in a separate 
PC. This PC also runs xPC Target software, but standalone 
configuration has be selected, because after design the model 
is rarely changed. 

Simulink model of the training ship Blue Lady is based on 
the set of non-linear equations, which describes dynamics of 
actuators dynamics, forces and moments acting the ship, hull 
dynamics and kinematics (Fig. 9) [1, 2]. Simulink model is 
created in a form of C-language S-function and also is used for 
computer of-line simulation of the first design stage (Fig. 3). 
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Fig. 5. Definition in Simulink of a decoder block 
of a NMEA 0183 format record from gyrcompass 


Simulink diagram prepared for Hardware-in-the-Loop 
trials (fig.10) includes some elements which fit parameters 
units from the Blue Lady model to those required by training 
ship communication interfaces (Fig. 2). Decoder and encoder 
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‘BL' in a received 
binary record from 
hub 


Verification of control 
byte in a binary 
record 
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Fig. 6. Definition in Simulink of a decoder block 
of a binary format record received from hub 


blocks in HiL diagram are revers counterparts of blocks in the 
sceleton. 

HiL stage drew up conditions for the simulation to the 
reality of the ship’s control. Circuit connections between the 
system control and object control are identical with those of the 
“Blue Lady”. In both cases the same Simulink block diagrams 
of the environment, the same compilers and starting procedure 
are used. Introduction of this stage saves a lot time for pre- 
testing of systems, debuging diagrams, verifing communication 
interfaces. 
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Fig. 7. Definition in Simulink of a encoder block for generation 
of binary format records transmitted to the hub 
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Fig. 9. Multidimensional Dynamics model 
of the training ship “Blue Lady” [1, 2] 
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Fig. 10. Simulator diagram of Blue Lady training Ship used for target xPC for Hardware-in-the-Loop trials 


CONCLUSIONS 


Matlab with its toolboxes Simulink, Simulink Coder, xPC 
Tatget is powerful software package for building a useful 
design environement in the area of a real-time control so Matlab 
enables posibility for rapid prototyping of control systems of 
ship movements. It is possible to build a multistage research 
(Fig. 3) test bed where proposed control algorithms can be 
tested in a computer simulation then during hardware-in-loop- 
simulations and finally during trials in a real environement. 
Very interesting option is replacing tests on real ships with trials 
in a reduced scale environement, when real ships are replaced 
by the reduced matrial ship models. The models are build in 
such way that rules of mechanical similarities are preserved 
so results for the reduced scale and real environement can be 
compared. 
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to Ship Trajectory Generation 


Andrzej Rak, M. Sc., 
Witold Gierusz, Ph. D., 
Gdynia Maritime University 


ABSTRACT 


This paper presents the application of the reinforcement learning algorithms 
to the task of autonomous determination of the ship trajectory during the 
in-harbour and harbour approaching manoeuvres. Authors used Markov 
decision processes formalism to build up the background of algorithm 
presentation. Two versions of RL algorithms were tested in the simulations: 
discrete (Q-learning) and continuous form (Least-Squares Policy Iteration). 
The results show that in both cases ship trajectory can be found. However 


discrete Q-learning algorithm suffered from many limitations (mainly curse 

of dimensionality) and practically is not applicable to the examined task. On the other hand, LSPI gave 

promising results. To be fully operational, proposed solution should be extended by taking into account 
ship heading and velocity and coupling with advanced multi-variable controller. 


Keywords: ship motion control; trajectory generation; autonomous navigation; 
reinforcement learning; least-squares policy iteration 


INTRODUCTION 


Development of control systems in recent years is 
determined by the expansion of relatively cheap computing 
power. This process is obviously present in the field of ship 
motion control systems too. As a result one can observe 
a progress from course keeping autopilots to trajectory tracking 
autopilots and further to integrated motion control systems with 
functionality of route planning, anti-collision subsystem and 
advanced multi-variable autopilot [3]. 
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TRACKING CONTROLLER 


Reference 
trajectory 


Desired 
position 


Current 
position 


ACTUATORS 


MEASUREMENT 
AND 


FILTERS 


Fig. 1. Block diagram of the proposed extension 
of the ship motion control system 


If the system is to be fully autonomous it should incorporate 
block of trajectory determination for in-harbour and harbour 
approaching manoeuvres. This task is significantly different 
from open waters or restricted areas anti-collision manoeuvres. 


The ship is moving with relatively small velocities and 
trajectory should be determined simultaneously for position, 
transversal and longitudinal velocities as well as heading 
signals. Trajectory generator has to cooperate with multi- 
variable autopilot capable to control all of mentioned signals 
[3, 4, 8]. The proposed location of such block in ship motion 
control system is shown in Fig. 1. 

Present paper describes the results of application of 
reinforcement learning (RL) algorithms to the generation of 
a ship trajectory. RL applications in the ship motion control 
system are quite rare. There are only a few examples in the 
openly accessible bibliographical sources [7, 9, 12]. 

In the first section the idea of the reinforcement learning 
and Markov decision processes (MDP) as a formal notation 
of RL problem are discussed, They deliver theoretical 
background of the algorithms. In the second and third sections 
learning algorithms in the discrete and continuous domains are 
presented. The results of computer simulations are shown in 
section fourth. Fifth section concludes the paper. 


REINFORCEMENT LEARNING AND 
MARKOV DECISION PROCESSES 


In reinforcement learning procedure a controller (agent) 
interacts with the process (environment) by means of three 
signals: a state signal which determines current state of the 
process, action signal which is used by the controller to 
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influence the process and a reward given by the reward function 
which gives the measure of the controller performance (see 
Fig. 2). In each consecutive time step the controller observes 
a state of a process and issues an action to move the process to 
a next state. At the same moment the reward function, based on 
a state and chosen action, evaluates this move issuing value of 
reward. In a next time instance whole procedure is repeated. 
The goal of the reinforcement learning process is to find (learn) 
a strategy of action selection for the controller which maximises 
the cumulative reward in a long term. 


Reward 
function 


state x 


Fig. 2. Interactions in the reinforcement learning process 


The dynamics of the process can be deterministic or 
stochastic. In this paper the deterministic case will be discussed: 
it means we assume that in particular state choice of certain 
action always gives the same next state). Extension to the 
stochastic case can be easy find in the references [1, 2, 10]. 

It is easy to notice that the problem of reinforcement 
learning can be described in the Markov decision process 
(MDP) formalism. Markov decision process is defined as 
a quadruple (X, U, f, p) where X = {x,, X» ..., Xpy is finite 
set of states; U = {u}, U, ..., Up} is finite set of actions and 
f: X x U — X is transitions function determining the state 
in a next time step: 


Xe = f(x, Uy) (1) 
At the same time the controller receives the value of reward 
according to the reward function p: X x U > [ R: 
Ik 7 Pio u,) (2) 
where we assume that ||p||,, = sup, „PC, u)| 


The controller chooses action according to its own policy 
h: X > U, using: 


u, = h(x,) (3) 


Taking into account above definitions we can state that 
given functions f and p as well as current state x, and action u, 
are sufficient to determine the next state x,,, and reward r,,,. 
This fulfils Markov property. 


LEARNING IN DISCRETE SPACE 


As we mentioned in the previous section, the goal of RL 
is to find an optimal policy that maximises the return from 
any initial state xp. The return is cumulative value of rewards 
collected along a trajectory originating in x,. There exist a few 
types of return definitions. We will use one of them; infinite 
horizon discounted sum: 


R”) = Hyp, hy] 4) 
k=0 


D This condition is partly violated by the need of exploration of 
the discrete RL algorithms [10]. 
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The infinite horizon return has better theoretical properties 
leading to the stationary optimal policies. But in a practical 
case the sum is limited by the number of steps in a trajectory 
between initial and final states. 

Discount factor y controls a trade-off between the quality of 
the solution and convergence rate of RL algorithm and usually 
is set by trial and error procedure. 

A convenient way to represent policies are their value 
functions [10]. In an area of RL two types value functions are 
used: state value function (V) and state-action value function 
(Q). The latter one is more general and incorporates the first. In 
this section algorithms based on the Q-functions are presented. 
It is a mapping Qh: X x U > [ R and represents a reward of 
choosing action u in a state x according to the followed policy 
h cumulated with the return from the next state [1]. 


h h 
Q' (x, u) = p(x, u) + yR [f (x, u)] (5) 
The optimal Q-function is defined as the best Q-function 
that can be obtained by any policy: 


Q* (x,u) = max Q" (x, u) © 


Any policy h* that selects at each state an action with the 
largest optimal Q-value i.e., that satisfies; 


h (x) earg maxQ (x,u) (7) 
h 


is optimal. In general, for a given Q-function Q, a policy h 
that satisfies: 


h(x)e arg maxQ (x,u) (8) 
h 


is said to be greedy in Q. So finding an optimal policy can be 
done by first finding Q*, and then using (7) to compute a greedy 
policy in it. Ifa process for which the learning is applied (Fig. 
2) is known (we have exact model of the process), Q-functions 
Qh and Q* can be easily found from the iterative Bellman 
equations [2, 10]. This leads to the dynamic programming 
Q-iteration algorithms [1]. 

In the case of RL the model of the process is unknown, so 
the next state and reward values are collected in the interactions 
with it. Therefore these algorithms are often called model-free. 
One of them, most widely used, is Q-learning: 


Qh & Ug) = Qh. (KU) + O) 


t EN +y maxQh Xk W)— Qk Xk Uk )] 


where a, € (0,1] is the learning rate. The term between square 

brackets is the temporal difference, i.e., the difference between 

the updated estimate r,,, + ymax Q,(X,,,,U’') of the optimal 
m 


Q-value of (x,, u,) and the current estimate Q(X% U,). 

As the number of transitions k approaches infinity 
Q-learning asymptotically converges to Q* if the state and 
action spaces are discrete and finite, and under the following 
conditions [11]: 

oO 


foe) 
- the sum Yo gives a finite value, whereas the sum Vo, 
0 ¥ k=0 


produces an infinite value. 
- all the state-action pairs are (asymptotically) visited 
infinitely often. 


The first condition is easy to satisfy. To fulfil the second 
one, a stochastic parameter ¢ is introduced. It represents the 
probability of selection of any action in encountered state. This 
is called exploration. On the contrary, the controller should 
also exploit current knowledge to improve performance by 
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selecting greedy actions in the current Q-function. The balance 
of exploration-exploitation is usually implemented in a form 
of e-greedy exploitation [10]. 


Algorithm 1. Q-learning with -greedy exploration 


Input: discount factor y, 
exploration schedule {g,} y 


learning rate schedule fag a n 


w 
ko? 
: initialise Q-function, e.g, Q} + 0 


— 


2: measure initial state xo 

3: for every time step k = 1, 2, ...do 
u € arg max,;Q,(x,,0) with probability 1- k 

4: Up & 


a uniformly randomin U with probability k 


5: apply ux, measure next state Xķ+ı and reward ry. 
Q, OK Uk) =Q Ok Uy) + 
+ oy [ey yt ymaxQ, (Xy pp Y) Z Q(X UI 


7: end for 


The complete algorithm, developed from equation (9) is 
presented in the frame Algorithm 1 in this section. This version 
of Q-learning was used in the computer simulations. 


LEARNING IN CONTINUOUS SPACES 
USING FUNCTION APPROXIMATORS 


Discrete algorithm for RL has a significant drawback. 
It requires an exact representation of a value function and 
policy. It means distinct values of the return estimates for 
each state-action pair has to be stored as well as actions for 
every state. When an action space or state space is large this 
can be extremely difficult or practically impossible. To reduce 
a number of parameters that has to be stored approximation 
techniques are used. 

Generally in the RL algorithms approximation is used 
not only for function representation. Policy iteration which 
will be introduced in the subsequent parts of this paper must 
repeatedly solve potentially difficult maximisation problems 
over the action variables (policy evaluation). This can be done 
by sample-based approximation. In this research parametric 
approximation was used. 

Parametric approximators are mappings from a parameter 
space into a space of functions [6]. The functional form of the 
mapping and the number of parameters are usually set by the 
skilled operator in advance and do not depend on data collected 
during the interaction with the process. The parameters are 
tuned using the data about target function. 

Let as consider the Q-function approximator parameterised 
by an n-dimensional vector 8. The approximator is the mapping 
F:R" — 3 where R” is the parameter space and 3 is the space 
of Q-functions. Every parameter vector 0 provides a compact 
representation of a corresponding approximate Q-function: 


Q(x,u) = [F@)](x,u) (10) 
where [F(6)](x,u) denotes the Q-function evaluated at the state- 
action pair (x,u). Therefore, instead of storing distinct Q-values 
for every pair (x,u), it is only necessary to store n parameters. 
But, it must be noticed that, since the set of Q-functions represent 
by F is only a subset of 3 any arbitrary Q-function, it can be 
reproduced only up to the certain approximation error [5]. 

The mapping F(@) can be generally non-linear. However, 
linearly parameterised approximators are preferred because they 
simplify an analysis of resulting RL algorithm. In presented 
algorithms we use linear approximators built with n Gaussian 


normalised radial basis functions (BF) 0,, ...0,: X x U>[R 
and n-dimensional vector of parameters 0. Approximate values 
were therefore computed as: 


FEJE Y= X o,&.ue,=6'& We an 
151 


Let us consider now model-based value iteration with 
parametric approximation. This procedure explains the ideas 
used in a final algorithm. One can observe that iteration formula 
(9) rewritten for model-based case can be generally stated as: 


Qr T T(Q,) (12) 


for consecutive | iterations, where T is Q-iteration mapping [1]. 
In approximate Q-iteration Q, cannot be represented exactly. 
Therefore, an approximation (10) has to be used: 


Q\= F@,) (13) 

Using this approximation in iteration formula (12) leads to: 
Qr = (Te FXO) (14) 

But the function T cannot be stored explicitly either. 


Instead it is represented by approximation using new parameter 
vector 0,,,. This vector is specified by the projection mapping 


P:5> R": 
0 PQ frn) (15) 


which should keep Q, 41 = F(0,,,) as close as possible to Qh, 
Usually a least-squares regression is chosen for P. Finally, 
approximate Q-iteration is a composition of mappings: 


0,.,=(P°T° F\6,) (16) 


The algorithm should be stopped when suitable parameter 
vector 0° is found. Eventually the estimate 6° should be kept 
as close as possible to the fixed point 6° of iteration (16). The 
whole procedure is illustrated on Fig. 3. The cycle of mappings: 
Q-value approximation (F), Q-iteration (T) and projection back 
to parameter space (P) is repeated until fixed point 0” is reached. 


(To F)(@,) 


T 
F(@,) 


parameter space 


Fig. 3. An idea of approximate Q-iteration 


Similar considerations lead to the analogous approximate 
policy evaluation algorithm for Q-functions. This algorithm 
starts from arbitrary chosen vector of parameters 05 and updates 
this vector in every iteration t using: 


0! = @ oT" oF)@") (17) 

As we mentioned before, approximators used in these 

considerations have linear properties. Based on this, it is 

possible to derive projected Bellman equation. Because state 

and action spaces are now finite by assumption we can rewrite 
policy evaluation mapping T”. 


[T"Q)(x, u)= 
p2 f(x, u, x’) [P(x, u, x’) + YOK, h(x’))] 


(18) 
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In a linear case approximate Q-function that has the form 
of (13) can be written as: 


Q'(x, u = o'(x, we" (19) 


where (x,u) = [0,(x,u), ..., ,(x,u)]? is the vector of BFs 
and 6° is the vector of parameters. This relationship satisfies 
approximate version of Bellman equation (projected Bellman 
equation): 


Ah hy (Ah 
Q= (PTQ) (20) 
where PY performs a weighted least-squares projection onto 


the space of approximate Q-functions spanned by the BFs. To 
derive proper algorithm let us rewrite (18) in a matrix form: 


T" (Q)= p+yfhQ 2) 
Using symbols @ to denote BF matrix and w to denote 
diagonal weight matrix of PY we can write: 


A 


Q= 98 (22) 
and Bellman equation (20) as: 
PYT"(Qt) = Q" (23) 
Rearranging this equation and substituting: 
[= 9'wo 
A=o'wfho (24) 
Z=Q' Wp 
the projected Bellman equation can be written in a final 
form: 
Tror =yA0 +z (25) 


Solving this equation leads directly to the Least-Squares 
Policy Iteration (LSPI) algorithm presented in a frame 
Algorithm 2 [6]. 


This procedure was employed in simulation software 
used to the generation of the ship trajectory in a continuous 
state space. The next section presents some of the simulation 
results. 


EXAMPLE OF SHIP TRAJECTORY 
GENERATION FOR NAVIGATION IN 
RESTRICTED WATERS 


Algorithm 2. Offline least-squares policy iteration 


Input: discount factor y, 
BFS ),...,,:XxU >R, samples {(x; ,u, .x' o dil, 51, n} 


initialise policy ho 


1: 

2: repeat at every iteration | = 0, 1, 2, ... 

3: Ty &0,Ay —0,z) <0 

4: forl,=1,...,n, do 

5 Di eD a HAE u 67 u) 
6: Ay SA a FOC u Mp" [xi hC )] 
7: Z & Z- HOX u, n, 

8: end for 


9: solve ir 0 = YLA, 0, +z. 
n, 3 n, So nn 3 


s s s 


10: h(x)—u, where u € argmax,o'(x,0)8);, Yx 
11: until h,,1 is satisfactory 
Output: h° =h,,, 
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The simulation experiments were done in the MATLAB 
environment. First the Q-learning algorithm was tested. Two 
variants were implemented: 

- Classic one: where the state transition can be done only 
to the neighbour-state in the discrete grid. It means that 
there are seven possible actions to chose. Six of them are 
“king moves” [10] to the neighbour state, and seventh is 
the loopback to the current state. The value function was 
set to -1 for all state transfers except last move to the goal 
state which yields 0 reward”). 

- Modified one: prepared to overcome limitations of 
neighbour-state only transitions. In this case action chosen 
can “jump” from every state to the any of states that are not 
blocked by obstacles. The reward is equal to the negative 
Euclidean distance from k + 1 state to the goal state. In this 
case action value function has to store n? distinct values, 
where n is the number of states in the whole grid. 


The results of the simulation for classic (step-by-step) 
variant are depicted on Fig. 4 and Fig. 5. Q-value function 
on Fig. 4 is marked by colours. More red for particular state 
means that the return (negative) associated with this state is 
bigger. It is easy to notice that trajectories from all states will 
tend towards the goal state; even starting point for the learning 
algorithm was not changed. This property is a result of e-greedy 
exploration part of the Algorithm 1. 


Starting point 


2 4 6 6 10 12 14 


Fig. 4. Q-value surface of Q-learning with step-by-step type strategy. 
Corresponding proposed trajectory is presented on the Fig. 5. 
Obstacles are dark blue. (2500 episodes) 


The corresponding greedy trajectory originating from the 
starting state (9, 14) is shown on Fig. 5. It can be easy noticed 
that it is not unique trajectory passing minimal number of state 
on the way to the goal. 

Results of the second variant (point-to-point) of the Q- 
learning algorithm are presented in Fig. 6. In this and next figure 
the colour map for values is reversed. Except modifications 
mentioned above, one more change to the standard algorithm 


2) Very often the reward is of negative value. It means it works as 
a “punishment”, not a “reward”. However it is a standard in RL not 
to change the term “reward” even for negative values. 
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was done. Along with the obstacle collision check in every step, 
there is a test if there exists straight, direct path from the current 
state to the goal. If so, then the state transition is done directly 
to the goal and algorithm finishes. This modification caused 
brown “shadows” over certain states of the grid. Because these 
states lie on the “open” side of the obstacles, learning algorithm 
did not “visit” them. The proposed modifications accelerated 
convergence of the solution. 


oe 
Pe ne See epee eee ae 
> è 
= te eer ee ane See eee oe 
r 
P A ORN AE anon E E ee ee 
ih ee cme were iene eee suspen 
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2 4 6 6 10 12 14 


Fig. 5. Final greedy trajectory of Q-learning with step-by-step type strategy 
- see Fig. 4 (2500 episodes) 


Starting point 


2 4 6 8 10 


Fig. 6. Results of Q-learning with point-to-point type strategy. 
Proposed trajectory is marked by the red doted line. 
Obstacles are white. (15 000 episodes) 


On the contrary to the step-by step variant this one produces 
unique, broken-line trajectory. Unfortunately, this version is 
exceptionally computing power demanding. On the Fig. 4 
there are 15 x 15 = 225 states in the grid. Multiplying them by 


7 possible actions yields 1575 state-action pairs for Q-function 
evaluation. The same grid for the second case will give 
(15 x 15)? = 50525 state action pairs to process. This number 
will expand very quickly if one wants to improve the precision 
of the ship positioning by the grid refinement or supplement 
a heading angle. 

It means that Q-learning algorithms in a discrete version 
are of limited usability for manoeuvring ship trajectory 
generation. 

Last figure presents the LSPI approximate Q-value surface 
(precisely: one of the cross sections through it) given by linear 
approximation over a (5x5) grid of RBFs. Additionally, action 
space was also approximated by (11x11) grid of RBFs. As 
one can notice, there is not the starting point for the learning 
algorithm. Version of the LSPI algorithm implemented in 
the research works in a batch mode doing calculations off- 
line for all previously collected data samples. The starting 
point was set to show only an exemplary trajectory. Whole 
surface is parameterised, therefore the problem of the grid 
scale mentioned earlier vanishes. The same factor makes it 
susceptible to heading control extension. 


Fig. 7. Least-Squares Policy Iteration. Surface of approximate 
Q-value function with proposed trajectory. 
Obstacle is white. (11 iterations for 17000 samples) 


CONCLUSIONS 


Results of the research presented in a previous section 
allow conclude that: 

- The discrete Q-learning algorithms can be used for ship 
trajectory generation only when precision of the motions 
is not of major importance or the maneuvring area is 
significantly small. 

- Point-to-point version of the algorithm is very vulnerable 
to the curse of dimensionality. 

- Off-line, continuous domain, LSPI algorithm seems to 
be good alternative to the discrete Q-learning. To test all 
its properties in a context of ship trajectory generation, 
it should be expanded to incorporate heading set-point 
values and, in anext step, longitudinal and transversal 
velocities. 
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Linear Matrix Inequalities in multivariable 
ship’s steering 


Monika Rybczak, M. Sc. 
Gdynia Maritime University 


ABSTRACT 


This paper explains the basics of the Linear Matrix Inequalities (LMI), with examples of 
simulations and calculations created in Matlab/Simulink programming environment where 
the controlled plant is the “Blue Lady” ship model. First chapter of this paper gives a short 
overview of publications describing the use of Linear Matrix Inequalities method. Second 
chapter contains basic definitions and equations for the LMI method. Chapter three presents 
the use of LMI method for ship control by describing controller synthesis for the “Blue 
Lady”. Chapter four compares the operation of two controllers, the first one consisting of 
three independent proper adjusted PID controllers and the second one being a multivariable LMI controller. 


Finally conclusions from the above comparison are drawn. 


Keywords: linear matrix inequality; use of LMI; synthesis of regulators; multivariable system; ship’s control 


APPLICATIONS OF THE LMI 


Linear Matrix Inequalities (LMI) methods are used mostly 
for theoretical numerical calculations. Additionally in the 
past few years arising demand for practical use of LMI in 
science and industrial applications, such as different kinds of 
complex networks (e.g. power systems, computer networks) 
and industrial engineering, was observed. 

The use of LMI for modelling synchronization problems 
in complex networks with rapidly changing topology was 
presented in [1]. Examples of such networks were the Internet 
and power grids. More specifically the LMI method was 
used to design a controller used for damping oscillations in 
a power grid which was described in [2]. The idea was to 
stop the controller operation when the grid was in a stable 
state. A controller with such a design could be used both in 
decentralized and centralized networks. Similar design was 
described in [3] where a controller was damping oscillations 
of asynchronous generator connected to a rigid power 
grid. A related example was also a thyristor controller for 
a capacitor bank stabilizing a flexible power grid which was 
described in [4]. Areal industrial application was described 
in [5] where a control system was controlling several steam 
boilers supplying a decentralized steam consumer network 
characterized with frequent and large load changes. A similar 
example was a controller for chemical to electrical energy 
conversion or more precisely a steam boiler — steam turbine 
with generator couple. The complexity of the process was 
caused by the need to parameter steam generation, steam 


input to the turbine and electrical energy from the generator 
parameters all at the same time as described in [6]. Another 
application was a system to control skyscraper movemen't 
under extreme wind conditions by the use of an active mass 
driver. The mass driver was installed on a real building model 
and was tested in a wind tunnel as described in [7]. 

The most closely related to marine applications, and thus the 
most practical use of LMI was the controller of an autonomous 
underwater vehicle MR-X1 (Marine Robot Experimental 1) 
built by the Japanese Agency for Marine-Earth Science and 
Technology (SJAMSTEC). The robot was equipped with 5 
propellers, it could work in remote control mode or as a fully 
autonomous unit. Aquarium test trials gave very promising 
results, a detailed description can be found in [8]. 


THE PRINCIPLES OF USE 
OF THE LMI METHOD 


Linear Matrix Inequalities LMI are described by canonical 
form [9, 10]: 


m 
F(x) =F, +% xF >0 (1) 
i=l 
where 
F(x) — variable affinic function of the variable x 
is a positive definite matrix, x E R™, 
F,e R™. 
F,,i=0, 1....m — constant and symmetrical, 
x — the variable (unknown). 
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The LMI method is related to the feasibility problem, which 
comes down to searching for the answer to the question whether 
there exists a solution x to the LMI problem in it’s overall form 
shown below: 


A(x) <0 (2) 


Where a is the state matrix of the control system. 

In order to create an LMI for a control system for the 
object is necessary to check if the eigenvalues of the matrix 
a of the controlled closed loop system are placed in the left 
half complex plane. 

Next the feasibility problem and stability can be checked 
with the Laypunov function shown below: 


V(x) = x! Px (3) 


Where P is a positively symmetric matrix. 
The stability condition can be formulated with the use of 
LMI as shown below: 


ATP + PA<0, P=P!>+0 (4) 


Fulfilling the P = PT > 0 (operand ,,>” means that 
eigenvalues of matrix P are positive) condition determines 
whether the control system is stable and solves the feasibility 
problem, which is finding a positive definite matrix P for the 
given state matrix A. After checking the above conditions 
dynamic properties of the control system can be designed 
by pole placement in a specific part of the left half complex 
plane. A defined plane for pole placement was designated 
C stab [12]. 

Three exemplary limit areas have been shown below and 
the first case has been used in simulations described in this 


paper: 
a) First case 


-r s+q 
C;tab ={S eC|s+q] <} <l- <0 (5) 
s+q -r 
Where: 
r — 
q = 


radius of the circle, 
centre of the circle 


Im 


Re 


Fig. 1. Left half plane with a circle 


b) Second case 


Cab = {S € Clo, <Re6)<a,}S 


= 6) 
s +5)-2a 0 ( 
( ) 2 _ <0 
0 - (s +s)+ 20) 
where: 
Q,,0, — vertical bars 


Fig. 2. Left half plane with two vertical bars 
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c) Third case 
Cyan = SE ClRe@)tg(0) < -imee 


G+ s)sin(0) Es +s)cos(0) 20 (7) 
(-s+s)cos(®) 6+s)sin(0) 


where: 
0<(0,5) — angle of flare 


Im 


Re 


Fig. 3. Left half plane with an angle of flare 0 sector 


After pole placement in the left half complex plane a further 
controller synthesis requires a description according to a defined 
standard such as H,,. 

A multidimensional control system MIMO (Multiple Input 
Multiple Output) for controlled system G is shown in Fig. 4, 


G 
Object 


Fig. 4. Structure of a system for tracking set value to synthesize 
a multivariable regulator 


Controller 


and has the below transfer function: 


Gy - 20. 20. 


(8) 
W(s) wœ 
where: 
Z(s) — Laplace transform of output vector, 
W(s) — Laplace transform of input vector, 
z(t) — output signal vector, 
w(t) — input signal vector. 


Controlled system G described by equation (8) was 
presented in matrix transfer operator form shown below: 


AB 
G= (9) 
C D 
For which the space state equalities are as follows: 
X= Ax + Bw 
(10) 
z = Cx + Dw 


where: 

A — state matrix with dimensions n x n, 

B  — input matrix with dimensions n xr, 

C — output matrix with dimensions m x n, 

D — direct feedback matrix with dimensions m x r, 

x — state of the system with a vector dimension n, 

w — input signal with a vector dimension r, 

z — output signal, with a vector dimension m, measured by 
sensors. 
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In control theory [9, 16] many criteria are used to determine 
the quality of control. In this paper the H,, standard was 
selected, it determines the relation between input signal vector 
w(t) and output signal vector z(t) for G controlled system 
(see Fig. 4): 


[GS], = SUP Sa lG) (11) 


Where o,,,,(Z) is the highest singular value of matrix A, 
which is equivalent to a square root of the highest eigenvalue 
of the matrix ATA. In numerical calculations the upper limit 
value of H,, was called scalar variable gamma y. 

The H,, standard for G system matrix is lesser than the 
scalar variable gamma y [12, 16] if, and only if y2I - DT > 0 
and there exists a matrix P= PT > 0 that fulfils the linear matrix 
inequality shown below: 


(ATP + PA + CTC) + (PB + CTD) 
(y2I — D™D)-1(B™P + DTC) < 0 


After calculations, the value of variable gamma is an 
approximation of the upper limit of the H, standard. During 
controller synthesis there are two possibilities, the H,, can 
be minimalized or can be a restriction if gamma is set to be 
a constant value. In this paper gamma, by means of trial and 
error, was set to be a constant value of 1.5. 

After using Schur’s complement to linear matrix inequality 
(12) we received the H,, standard in the form of an LMI 
requirement: 


(12) 


A'P+PA+C'C PB+C'D 
aai 0 (13) 
B'P+D'C D'D-y’I 
If we assume that the transmittance T,___.(s) describes the 


control system deviation transmittance, where z,, is defined as 
Z„ = W —z (see Fig. 4) then the minimalization of the standard 
of the defined transmittance tends to assume the smallest value 
for all frequencies. 


CASE STUDY 
a) Controlled system description 


The training ship “Blue Lady” is a 1:24 scale model of 
a VLCC (Very Large Crude Carrier) tanker that was meant to 
carry large amounts of crude oil but unfortunately was never 
built. 

Main propulsion is a d.c. electric motor driving a five 
blade, fixed pitch, propeller through a reduction gear [15]. The 
propulsion characteristic corresponds to turbine propulsion but 
it has an option to be changed to diesel engine propulsion. The 
model is also equipped with two tunnel thrusters and two side 
thrusters simulating 60 ton tug boats. ,,Blue Lady’s” cockpit 
has enough space for two people. The silhouette of the ship is 
shown below. 


GPS antenna 


Fig. 5. Silhouette of “Blue Lady” with cockpit arrangement 
and GPS antenna 


Tab. 1. “Blue Lady” model and real ship parameter comparison 


Overall lenght 
Breadth 


330.65 [m] 
57.00 [m] 
28150[m/h] 


13.78 [m] 
2.38 [m] 
5741.2 [m/h] 


Maximum speed 


Training ship “Blue Lady”, being the controlled system 
in the considered case, had three input signals [t, Ty Ty] (see 
Fig. 6) where: 


T, — required force (thrust) in the ships longitudinal axis, 
T, — required force (thrust) in the ships lateral axis, 
T, — required turning moment. 


Taking into consideration the number and type of propellers 
eight command signals for propulsion and steering equipment 
were implemented [ngc dc sstd, sstr, ssod, a4, ssor, a,,] (see 
Fig. 6). 

Where: 

[X Y, N,] — were forces and moments created by propulsion 
and steering equipment, and the three output 
signals were position coordinates x(t), y(t) and 
the heading y(t). 


For future controller synthesis with the LMI method the 
state model aspect had to be determined. Considering input 
and output signals of the future controller for the “Blue Lady” 
model the following were taken into account: dynamics and 
kinematics of the ship, Kalman’s filter, geographic coordinates 
and a system for recalculating command signals to particular 
thrusters. It should be noted that dynamic equations of the ship 
had the following form (so called 3DOF model): 


m(ù- rv) = X,[N] 


m(v—ru)= Y,[N] (14) 
Izr =N,[N] 
Where: 
m — ship’s displacement, 
u — linear velocity in respect to X axis, 
v — linear velocity in respect to Y axis, 


r — angular velocity around Z axis. 

Whereas course changes and ships position which were 
the kinematic equations had been calculated from the below 
formula: 

t 


W(t)= Volto) + rode 


t 
t 0 


X(t) = xolt) + [u cos(y()) -vsin odr 


to 
t 


y(t) =yo(to) + jo sin(\y(t)) + v(t) cos(y(t))]dt 
to 
Ablock diagram of the controlled system has been presented 
on Fig. 6. 


(15) 


b) Linear model of the controlled object 
It turned out, during identification process, that three signal 


channels demonstrated weak correlation between output and 
input signals: 
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Fig. 6. Block diagram of the controlled system for identification process. 
u, v, r — ships velocities; x, y, y — ships position and course; 
[ng., 6... ]T vectors — command signals for propulsion 
and steering equipment, [X, Y; N,]T vectors — forces 
and moments generated by propulsion and steering equipment 


N N A (16) 
T> V,Ty > UT, >I 


therefore these subsystems were canceled from the final model 
(see Fig. 7): 


Fig. 7. Control channels after identification process 


Finally the multivariable linear state model of the system 
had the following form: 


X,| Jau 0 0 X, 

X,/=| 0 ay, a,/*|x.|+ 

X 3 a ur a yr a ir, X3 
bw 0 0 Tg 

+| 0 by ba” Ty (17) 

bir byr De Ty 
u 1 0 0 Xij 
v|=|0 1 O]*I|x, 


r 0 0 1 X3 


where 0 values denoted cancelled channels. 


The coefficient values of the model parameters were 
obtained as average values from all identification experiments. 
State model coefficients were as follows [11]: 


Tab. 2. Values of ”Blue Lady” model coefficients in state space 


parameter 


parameter 
3.62 * 103 
2.06 * 10° 
1.61 * 10° 
3.00 * 10° 
1.15 * 10° 
7,00 * 103 


-3.36 * 103 
-9.00 * 103 
-2.00 * 104 
-3.00 * 103 
-1.00 * 103 
-7.75 * 10° 
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c) Controller synthesis 


The first condition for the LMI method applied to the 
feasibility problem (2) and required checking ifthe eigenvalues 
of the system matrix were placed in the left half complex 
plane. For the state matrix of the simulated system a that had 
the below form: 


aw 0 0 
A=| 0 awy ây (18) 
aur ay, air 


Next the stability had to be checked, which came down 
to fulfilling the Laypunov inequality (4) which meant that for 
the given state matrix a you looked for a positively determined 
symmetric matrix P = PT > 0 [9, 12, 13]. 

For the state matrix of the simulated controller a symmetric 
matrix P was found, and had the following values: 


2.0658 0.0119 0.1021 
P=/0.0119 1.9766 —-—0.0044 (19) 
0.1021 -—0.0044 1.9034 


After describing LMI conditions from the theoretical part 
of the article with matrix form, a multivariable controller was 
designed for a MIMO system, the described ”Blue Lady” 
ship model. First step was to determine the matrices based 
on system structure (see Fig. 4) and system description (22). 
Below a system of equations for the controlled system has been 
presented together with output equations defining signals for 
H,, standard, measured signal from system to controller “e” 
and output signal “z”. 


A B By 
Co Dou Dow u 
zZ |= $ (20) 
C, Des Dw Ww 
e 
Ce Day Dew 


LMI matrices had the following form: 


X=Ax+ B u +B, w 


Za = CX + Dout +D oy W 21) 
z=C,x+D,,u+D,,w 


e = C.X + D au +D ey W 


where the matrices representing input signals ”w” and ”u” 
were vector matrices [3x1] because they applied to [t, , Ty » Tpl 
signals. Additionally vector matrices for output signals “z” 
and ”z,,” were also [3x1] and they applied to: longitudinal “u”, 
lateral “v” and rotational “r” velocity signals. 

After checking that the controlled system fulfilled the 
feasibility condition (2) and stability condition (4) its dynamic 
properties had to be specified by pole placement in the left 
half complex plane [12, 14]. In this paper limit area (5) was 
used. 

For the input and calculation of the above parameters 
with the LMI method the LMI Control Toolbox in Matlab 
together with additional toolboxes SeDuMi (Self — Dual 
— Minimization) and YALMIP (Yet Another LMI Preprocessor) 
was used [18, 19]. 
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Fig. 8. Up side is the structure of a system for tracking set value, down side 
is the standard structure using the LMI method to synthesize a multivariable 
controller where the object is “Blue Lady” ship model 


Experimentally selected parameters of the limit area had 
the following values (see Fig. 1): 

- radius of the circle — 0.04 

- centre of the circle — 0.02 


For controller synthesis, after determining the limit area, 
additional evaluation of the H,, standard was done, for which 
the scalar variable “y” was the estimated upper limit value. 

After using Schur’s complement and making some 
modifications to the controlled system described by (24) and 
based on (13) the following H,, standard condition for LMI 
method was received: 


AR+RA'+B,C,+C/B, AL +A+B,D,C, 
A,+AT+C/DIB! A'S+SA+B,C, +CIB} 
Bi +D1 DIB! BIS+D B! 
CR +Da C, C, +Da D;C, = 
B,+B,D,D RC, +C Da 
SB,+B,D,, G" +CiDID!, 
| D+ DD 
Dow + Dau DDow -yI 


In the above matrix inequality new matrices R and S were 
used. They were calculated by congruence transformation of 
LMI conditions which was true for the below equations of 
a closed loop circuit: 


D, =D, 
C, =C,M" +D,C,R 

E (23) 
B, = NB, +SB,D, 


A, =NA,M" +NB,C,R +SB,C,M" +S(A+B,D,C,)R 


Above equations contained N and M matrices which were 
connected to Laypunov variable factorization of matrix X 
related to LMI conditions describing H, standard together with 
pole placement limit area, as described in [9, 12, 16]. 

For the simulations in this paper H,, was calculated to be: 
1.08, for a constant upper limit value of the scalar variable “y” 
equal to: 1.5. After calculating LMI method conditions the 
matrices for the designed controller were following: 


aral 0.0001 0.0016 -0.0196 0.0291 - 0.0041 - 0.0007 


-0.0276 0.0019 -0.0001 0.0031 0.0178 0.0708 | 
0.0024 - 0.0301 0.0015 - 0.0072 -0.0424 0.0098 


0.0007 0.0010 - 0.0058 - 0.0277 - 0.0035 0.0005 
0.0030 0.0077 0.0014 - 0.0031 -0.0426 0.0014 
- 0.0014 - 0.0047 0.0001 0.0007 0.0022 - 0.0490) 


- 0.0143 0.0000 0.0010 
0.0045 -0.0014 0.0153 
- 0.0003 - 0.0226 - 0.0010 
0.0040 0.0623 0.0225 
0.0211 -0.0181 0.0999 
0.0855 0.0020 - 0.0517 


(24) 


- 0.3026 0.0489 - 0.0035 0.0552 0.3187 1.1398 
0.0003 -0.0141 -0.2773 0.8015 -0.2311 0.0231 
0.0193 0.3300 - 0.0301 0.1415 0.6626 - 0.4161 


Ck 


0 0 0 
Dk=/0 0 0 
0 0 0 


SIMULATION RESULTS 


For the simulations in this paper basic assumptions have 
been analyzed. Fig. 9 below shows changes made to the 
structure of the system from Fig. 8 


e u - / 
Controller G 
7 - i el Object | 


Object 


Controller 


Fig. 9. Upper is the structure of a system for tracking set value with 
dynamic and kinematic parameters, lower is the standard structure using 
the LMI method to synthesize a multivariable regulator where the system is 
“Blue Lady” ship model 


The practical aim was to realize ”Blue Lady” ship 
model steering for three velocities. In the simulations 
a multidimensional LMI controller designed by the author was 
compared with a controller consisting of three independent 
PID controllers designed by an M.Sc. graduate of Maritime 
Academy in Gdynia, as described in his paper [17]. Simulation 
results were achieved in Matlab with YALMIP and SeDuMi 
toolboxes and had parameters from the table 3. 
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Tab. 3. simulation parameters 


Simulations with one velocity 
longitudinal | lateral 
velocity velocity 
[m/s] [m/s] 


rotational 
velocity 
[rad/s] 


simulation no. 


Simulations with two velocitie 
lateral 
velocity 


rotational 
velocity 
[rad/s] 


longitudinal 
velocity 
[m/s] 
0.2 
-0.210.2 | -0.08 10.08 
Simulations with three velocities 
longitudinal | lateral 
velocity velocity 
[m/s] [m/s] 
0.08 
-0.08 1 0.08 


simulation no. 


rotational 
velocity 
[rad/s] 


simulation no. 


Every simulation was started from yọ = 0 [deg]. 
In order to determine the quality of control the following 
formula was used: 


1 n 
Ax? =>) (X-X) (25) 
ni 
where: 
n — number of measurements, 
x — reference signal value, 
X,.— Output signal value received from the system. 


a) Simulations for one velocity, with a given value of 
longitudinal velocity u = 0.2 [m/s] 
Symulation with MIMO controler LMI 


o eyo 
0 E 0 


Symulation with controller PID 
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Fig. 10. Graphs for three components of sideways movement. Solid line 
— received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results 


Tab. 4. Comparison of average value deviations from given velocities 


A for PID for LMI 
0.0004433 0.0009162 
0.0000000 0.0000010 


0.0000004 0.0000501 


longitudinal — u [m/s] 
lateral — v[m/s] 
rotational — r [rad/s] 
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Comment: 

For the LMI controller the average deviation value 
from given longitudinal velocity was twice as big as for the 
PID controller. Settling time for the received velocity was 
comparable for both controllers. Unfortunately the average 
deviation value from given rotational velocity was significantly 
higher for the LMI controller than for the PID controller. 


b) Simulations for two velocities, with given values of 
longitudinal velocity u = 0.2 [m/s] and lateral velocity 
v = 0.08 [m/s] 


Symulation with controller PID Symulation with MIMO controler LMI 
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Fig. 11. Graphs for three components of sideways movement. Solid line 
— received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results. 


Tab. 5. Comparison of average value deviations from given velocities 


for PID 
controller 


for LMI 
controller 
longitudinal — u [m/s] 

lateral — v[m/s] 
rotational — r [rad/s] 


Comment: 

Average deviation value from given lateral velocity for 
the PID controller was unacceptable since the settling time 
was significantly higher than 2000 [s]. The deviation factor 
for lateral velocity was more than 30 times smaller for the 
LMI controller than for the PID controller. The remaining 
deviations for both controllers were considered standard and 
acceptable. 


c) Simulations for three velocities, with given values 
of longitudinal velocity u = 0.2 [m/s], lateral velocity 
v = 0.08 [m/s] and rotational velocity r = 0.1 [rad/s] 


Tab. 6. Comparison of average value deviations from given velocities 


for PID 
controller 
0.0004 
0.0041 
0.0001 


for LMI 
controller 
0.0010 
0.0002 
0.0002 


Velocity 


longitudinal — u [m/s] 
lateral — v[m/s] 
rotational — r [rad/s] 


Comment: 

Comparing both controllers for this case the LMI controller 
had better parameters. Lateral velocity settling time was 
significantly shorter (approx. 1000 [s]) compared to PID 
controller (more than 2000 [s]). Deviation factor for lateral 
velocity was almost 40 times smaller for the LMI controller than 
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for the PID controller. For longitudinal and rotational velocities 
both controllers had comparable parameters. 
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Fig. 12. Graphs for three components of sideways movement. Solid line 
— received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results 


d) Simulations for one velocity, with a given value of 
longitudinal velocity alternating between u = 0.2 [m/s] 
and u = - 0.2 [m/s] in time instant equal to 0[s], 600 [s], 
1100 [s] and 1600 [s] 


Symulation with controller PID Symulation with MIMO controler LMI 
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Fig. 13. Graphs for three components of sideways movement. Solid line 
— received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results. 


Tab. 7. Comparison of average value deviations from given velocities 


for PID for LMI 


controller controller 
longitudinal — u [m/s] 0.0124 0.0152 


lateral — v[m/s] 0.0000 0.0000 
rotational — r [rad/s] 0.0001 0.0019 


Comment: 

Above simulations presented both controllers with given 
longitudinal velocity changes during the simulation. Average 
deviation value from given longitudinal velocity was similar for 
both controllers. The LMI controller however had worse results 
for given zero velocities, there have been some interference in 
the system. 


Velocity 


e) Simulations for two velocities, with given values of 
longitudinal velocity alternating between u = 0.2 [m/s] 
and u = -0.2 [m/s] and lateral velocity alternating 
between v = 0.08 [m/s] and v =- 0.08 [m/s] in time instant 
equal to 0 [s], 600 [s], 1100 [s] and 1600 [s] 


Symulation with controller PID Symulation with MIMO controler LMI 
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Fig. 14. Graphs for three components of sideways movement. Solid line 
— received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results 


Tab. 8. Comparison of average value deviations from given velocities 


for PID 
controller 

0.0127 

0.0041 


for LMI 
controller 
0.0133 
0.0030 


Velocity 


longitudinal — u [m/s] 
lateral — v[m/s] 


rotational — r [rad/s] 0.0001 0.0024 


Comment: 

In the above case both controllers were presented with given 
longitudinal and lateral velocity changes during the simulation. 
Looking at settling times and average deviation values from 
given velocities the LMI controller had worse results than in 
the previous case but they were still satisfactory as for this 
stage of the design. 


f) Simulations for three velocities, with given values of 
longitudinal velocity alternating between u = 0.2 [m/s] 
and u = - 0.2 [m/s], lateral velocity alternating between 
v= 0.08 [m/s] and v =- 0.08 [m/s] and rotational velocity 
alternating between r= 0.2 [rad/s] and r = - 0.2[rad/s] in 
time instant equal to 0 [s], 600 [s], 1100 [s] and 1600 [s] 


Symulation with controller PID Symulation with MIMO controler LMI 
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Fig. 15. Graphs for three components of sideways movement. Solid line 
—received velocity, dashed line — given velocity. Left side — PID controller, 
right side LMI controller simulation results 
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Tab. 9. Comparison of average value deviations from given velocities 


for PID 
controller 
0.0125 
0.0040 
0.0011 


for LMI 
controller 
0.0153 
0.0031 
0.0076 


Velocity 


longitudinal — u [m/s] 
lateral — v[m/s] 
rotational — r [rad/s] 


Comment: 

In the above simulations both controllers were presented 
with given longitudinal, lateral and rotational velocity changes 
during the simulation. Best results for the LMI controller were 
obtained for deviation values from given lateral velocity, for 
the other two velocities the PID controller presented better 
results. The LMI controller however, had better settling times 
than the PID controller. In the LMI controller, despite having 
a worse result for rotational velocity, the settling time was not 
affected as much as the settling time for lateral velocity in the 
PID controller. 


CONCLUSIONS 


This paper was based on results of simulations for a “Blue 
Lady” ship model trajectory simulator created in Matlab/ 
Simulink. For LMI calculations to be possible additional 
Yalmip and SeDuMi toolboxes had to be installed into Matlab 
libraries. Above simulations have proven that both the LMI 
and PID controllers have different proprieties. In order for the 
correct controller to be selected the brief foredesigns have to 
be taken into consideration. If the important factor would be 
the output signal settling time than the LMI controller would 
be the better choice. Despite having different average value 
deviations from given signals the LMI controller was more 
resistant to interference during the simulations. Furthermore 
the LMI controller was of the sixth stage and the control was 
multidimensional. In future studies external interference will be 
added to the system to allow further study of the controller. 
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The branch-and-bound method and genetic 
algorithm in avoidance of ships collisions in 
fuzzy environment 


Mostefa Mohamed-Seghir, Ph. D. 
Gdynia Maritime University 


ABSTRACT 


Marine navigation consists in continuous observation of the situation at sea, determination 

the anti-collision manoeuvre. So it necessary to determine ship safe trajectory as a sequence 

of ship course changing manoeuvres. Each manoeuvre is undertaken on the basis of 

information obtained from the anti-collision system ARPA. This paper describes a method 

of safe ship control in the collision situation in a fuzzy environment based on a branch 

and bound method and a genetic algorithm. The optimal safe ship trajectory in a collision 
situation is presented as multistage decision-making process. 


Keywords: safe ship control; optimal control; safe trajectories; a branch and bound method; 


genetic algorithm; ship control; fuzzy set theory 


INTRODUCTION 


The research for effective methods to avoidance ship 
collisions has become important with the increasing size, speed 
and number of ships participating in sea transport. Since when 
were debut the application of the ARPA the safety of maritime 
navigation has increased. 

A new tendency within the contemporary domain of 
ship control involves an automation process of selecting an 
optimal manoeuvre or optimal safety trajectory, based on the 
information obtained from the anti-collision system [3, 4, 5, 
8, 9]. This paper discusses in detail the ship optimal position 
determining process, involving stages of ship trajectory, 
based on the kinematic model. It is assumed that the motion 
of targets is uniform and occurs as a straight line. Due to the 
fuzziness of the process, individual approach of a particular 
officer-navigator, the decision-making process is also, to some 
extend, an ambiguous evaluation of the safe approach distance 
and safe time to avoid collision maneuver. Moreover, it is 


Unfuzzy 
information 
from ARPA 


Fuzzyfication 


Set of 
principles 


assumed that an optimal safe trajectory in a collision situation is 
a multistage decision-making process in a fuzzy environment. 
The maneuverability parameters of the ship as well as the 
navigator’s subjective assessment are taken into consideration 
in the process model figure 1. 


PROCESSES OF SAFE SHIP CONTROL IN 
FUZZY ENVIRONMENT 


In order to describe the safe ship trajectory, the movement 
of a ship returning by her rudder in deep water was described. 
Still, the safe ship trajectory may prove little inefficient to 
evaluate the ship properties. Therefore in order to precisely 
calculate a ship dynamic properties we used parameters of the 
transfer function or the advance time t, and maximum angular 
speed  [2, 5, 6]. 

The model of safe ship trajectory can be represented by 
the state equation: 


f(X,8) > XxS (1) 


Unfuzzy result 

for supporting 

navigator’s 
decision 


Defuzzyfication 


Fig. 1. Block diagram of application of fuzzy sets in a safe ship control 
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Xa AFRL An (2) 


where: 

Xi XE X= {ap ap, ap-1, py 3p at — set of real 
ship position 
co-ordinates 

S E U= {Cos Cis -> Cat — control set 

The process comes to an end when a ship attains back points 

(final points) called the final states W C X [1, 2]. 


W = {ap ap+2? a} (3) 
The set of final states must satisfy this condition: 
Copt UR = LR safe (4) 
where: 
Copt (Wop Vopt) — Optimal control, 
Lp — membership function of fuzzy set collision 
risk. 


This membership function of fuzzy set collision risk can 
be presented in the form Figure 2, [2, 3]: 


ZCoxXxxX (5) 
Hg: XxX—[0,1]JER (6) 
: l 
Lr (k,j) = (7) 


EXP(Agis(K,j)CPA + Aer (Kj) TCPA}) 


HR 


CPA, 


Fig. 2. Graphical interpretation membership function 
of fuzzy set collision risk 


Similarly, the membership function of fuzzy set of goal can 
be written in the form: 


Gcxxs (8) 
gi XxU—[0,1JeR - 

l 1 
ee 10 


Now, it the membership function of fuzzy set constraints 
must be defined as constraints of maneuver at each step: 


Ccxxs (11) 
uc X x U > [0,1] € R (12) 
: (13) 


Belk) = pi AW) Veosy(k) — Veosy(k— DIE 
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The fuzzy set decision is determined as the fuzzy set 
DcXxsS. It isa result of an operation “*” of the fuzzy set of 
a goal and fuzzy set of constraints: 


D=G+*C (14) 
Hp (>) = He Ga) * Ha (>) (15) 

where: 

hap Agp Ag Ae — navigator’s subjective parameters, 

CP. — the Closest Point of Approach, 

TCPA — the Time to Closest Point of Approach, 

V — ship speed, 

Y — ship course. 


In order to find a solution to this the author proposed the 
following methods. 


BRANCH-AND-BOUND METHOD 


A fuzzy decision is a result ofa certain compromise between 
these sets G (fuzzy set of goals) and C (fuzzy set of constraints), 
if the trajectory is called sequence states attained, then 
membership function of fuzzy set decision define as [3, 4]. 


UplSo Sys +2 Sal Xo = 
= WelSo) * WgX,) * He(Si) Hg (X) * 16) 
Ti Uc(S,_1) i Ug (X,,) 


where: 


A — minimum. 


To maximise a membership function of fuzzy set decision, 
at using minimum type, we obtain the optimal decision [1, 5] 


Hp(Sp°) S1", e Sy.) Xo) = 
= max (iMe(SlXo) A Hg (SIX) 0D 
ae ae 


n- 


Kado eS. 4) ^ AX 


The decision process can be conveniently represented in the 
form of a decision tree, the root of the tree is the initial state 
Xo- We start from xo, and looking for the optimal decision, after 
that to put it to control uy and we pass to state x,. We determine 
again the optimal decision u, and we pass to the next state, until 
we attain the final state. In this manner we obtain sequence of 
states which present a ship’s optimal safe trajectory [4]. 


No = He (So) A UG (X,) 


Ny = Ue (So) A UG (Xo) A Uc (S,) A Ug (X) z 
=No A [uc (S,) A UG (X,)] 


Nk = [Me (So) AUG (X,)] A [He (S,) AUG (X,)] A 
A [Me (Sy) A Bg Ken) = (18) 
=k A [uc (Sp ^ Ug Xen] 


Np-1 ~ [uc (So) ^ Ug (X,)] A [uc (S,) AUG (X,)] A 
A [uc (Sa) A UG Xp] = 
=o A [He (S11) A Wg Xp] 


If range controls a mount S}, S,, ..., S,, K < n-1, it for each 
L>k, L<n-l it gets 


(19) 


1h 2 Ny 
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From this inequality emerge, that each value n; cannot be 
greater than value n, in other case at use of operation minimum 
A”. This way, it is possible to ascertain progressing, that 
inequality gets: 


Nk Z We = HplSo Si -> Sega 20 


To suppose, that it obtain k control stage and certain state of 
process state, now it must be choose optimal state from states 
achieved earlier. 

To continue this procedure until we obtain final state, 
the process ends and we get optimal safe ship’s trajectory in 
collision situation. 


Beginning of 
calculation? 


Loading 2 from the file for the 
selected ship size and the specified 
visibility 


Determination and memorizing the 
distance in which the ship is no longer 
endangered 


Calculation of the trajectory o 
foreign objects at particular 
stages 


Creation of a random way initial 
population of springs 
controls P(0 


Determination of the condition 
X of the ship at particular stages of 
their populations controls P(0) 


Evaluation of the initial 
population P(0) 


Sorting of 
population by descending 
according to the evaluatio 
unction 


Reporting the 
calculation 
method 


Initializing of variables, 
arrays and structures 


Crossbreeding of the T(t) and 
Creation of the population O(t) 
Mutation of the population of O(t) 


GENETIC ALGORITHM METHOD (AG) 


In the context of the multistage fuzzy control, an individual 
approach is understood as a sequence of particular control 
stages Sp, ...,S,,,. According to this approach of an individual, 
the fuzzy environment is assessed using a membership function 
of fuzzy decision. On the basis of this decision several potential 
solutions can be selected. The set of these solutions is termed 
as the population. It is assumed that the algorithm will operate 
on a population of a certain size, which is initially randomly 
generated. Some members of the population, which play the 
role of “parents”, are reproduced by crossover and mutation, 


(t<max number 
of generations) and 

(is improvement) and 
(uD<given threshold) 


Ranked or proportional reproduction of 
P(t) and creation of the population T(t) 


Succession — an elite selection 


P(t+1) := Pt) + OV(t) 


Determination of the condition X of the ship 
at particular stages of their population 
controls P(t +1) 


Evaluation of the new, 
base population P(t +1) 


population by descending 
according to the evaluation 


Fig. 3. The algorithm to determine the safe ship trajectory using AG 
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whereas new solutions are ,,children’”, The best and strongest 
of them all “survives”, this means that it participate further in 
this process. At the end of the process one can expect to find 
a very good, or perhaps even the optimal solution. 

Determining a safe trajectory of a ship, as formulated 
above, and finding the optimal sequence can be done through 
[1, 2, 6, 7]: 


Hp% Si» : Ae) = max Lp(Sp.S), . »Sn1|Xo) = 


ae ae 
ESD AXD A WAS) A WG (XQ) A... D 


max n— n 
SoS eA BES, Da AX 


Where at each stage (t-1), control S,, € U, the fuzzy 
constraints 1'5'(S,_,) are imposed, and at stage t on state X, the 
fuzzy goals ut (X;) are imposed. 

Before a genetic algorithm can be used some assumptions 
must be made as follows: 

- the problem isa series of controls So, ..., Sy} due to the relief 
of genetic algorithm and simplification of the results analysis. 
We assumed that a ship speed is constant and the control 
S,,, at stage t +1 is defined as the angle of the course Y, 
relative to the previous angle of the course Y, at stage t, 

- the actual encoding does not change. We assumed that 
the real representation of each gene S, € [0.360] is evenly 
distributed, which is natural in this case, 

- purpose function of each individual is a membership 
function of the fuzzy decision type, of minimum (A) given 
by the formula, 


Lp (So.8,. 79 Saa | Xo) = 
= ue (So) auo (X aue (S) A HE (Xp) A... 
a. ABE (Sn) ALG (Xn) 


- individuals with the highest value of the purpose function 
have the largest share in the next parental population, and 
weaker individuals are rejected in the selection process, 

- selection process is reproduction - proportional or 
ranked, 

- cross-sectional averaging is used as the most suitable for 
encoding floating point, 

- there occur perturbations of gene mutations, acting in 
accordance with the Cauchy distribution with a certain 


(22) 


a) 15; — 
Y[nm]| | Optimal trajectory 


105 


-5 + 


-15 -10 -5 0 5 X[nm] 15 


probability of occurrence in successive generations, 

- creating the next generation does not reject all parental 
individuals of the population, but the best of them attach 
to the descendants of individuals, 

- stopping conditions should involve at least the assumption 
that the algorithm terminates the calculation when the 
improvement is smaller than a certain threshold. 


SIMULATION RESEARCH 


The exemplary results of the performed computer 
simulations are as follows. 

The first, the computer simulation involved 20 targets. 
Through the adjustments we can determine whether the 
designated safe trajectories are indeed optimal (Fig.4). 


Tab. 1. 
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Fig. 4. The result of simulation of the collision situation in passing with 20 moving targets. a) Branch-and-Bound Method, b) Genetic Algorithm 


48 POLISH MARITIME RESEARCH, Special Issue 2012/S1 


The branch-and-bound method and genetic algorithm in avoidance of ships collisions in fuzzy environment 


In the case of the method of branch-and-bound. The 
ship holds his course given only after 24 minutes made the 
manoeuvre rate (y = 345°) to the left, after this stage, returns to 
the set course. After 6 minutes, the ship performs again a few 
manoeuvres to the left avoiding the last conflict of the object 
(9) before returning to a given course. 

In the case of the method of branch-and-bound. 

In this situation, the manoeuvre was performed on the right 
course (y, = 12.5°), the own ship returns to a given course 
by making some quick manoeuvres to the left passing few 
targets. 


CONCLUSIONS 


- In this paper it presented tow methods based on fuzzy 
set theory to solve task of optimal safe ship trajectory in 
collision situation, according to the international Rules of 
the Road at Sea. 

- This work showed a possibility of using a branch-and-bound 
method and genetic algorithms in a fuzzy environment as 
a method to solve the problem related to determining an 
optimal safe ship trajectory in collision situations with the 
use of computer software. 

- Itis possible to effectively solve tasks of determining a safe 
trajectory ofa ship as a multistage decision-making process 
in a fuzzy environment. 

- A navigator controls the ship with according to his 
individual assessment of the risk collision. This is an 
individual decision, which is not only subject to a condition 
of passing a greater distance from its set point. 

- This present paper showed that the suggested idea of 
a fuzzy set theory application is a promising way to solve 
the considered task and design novel anti-collision systems 
in the future. The fuzzy set theory can be applied in many 
domains. 


BIBLIOGRAPHY 


1. Kacprzyk J.: Wieloetapowe sterowanie rozmyte. WNT, 
Warszawa 2001 

2. Lisowski J., Mohamed-Seghir M.: The Safe Ship Control with 
Minimum Risk of Collision. Wit Press Comp. Mech. Publ., 
Boston 1998. 

3. Lisowski J.: Mathematical modeling of a safe ship optima 
control process. Polish Journal of Environmental studies, 14 (D, 
68, 2005. 

4. Mohamed-Seghir M.: Wybrane metody rozwiązania zagadnienia 
bezpiecznego sterowania statkiem w sytuacjach kolizyjnych w 
rozmytym otoczeniu. Wydawnictwo Akademii Morskiej, Gdynia, 
2007. 

5. Mohamed-Seghir M.: Optymalna bezpieczna Trajektoria statku 
w Rozmytym Otoczeniu. XIII Międzynarodowa Konferencja 
Naukowo-Techniczna, Gdynia 2002. 

6. Mohamed-Seghir M.: Methods based on fuzzy sets to solve 
problems of Safe Ship control. Novel Algorithms and Techniques 
in Telecommunications and Networking, USA Springer 2010, p. 
373-377. 

7. Mohamed-Seghir M.: The safe ships control in fuzzy 
environment using a genetic algorithm. Solid State Phenomena. 
Vol. 180, Trans Tech Publications, Switzerland, 2012, p 70-75. 

8. Pietrzykowski Z.: The navigational decision support system on 
a sea-going vessel. Maritime University, Szczecin, 2011. 

9. Sztapezynski R., Smierzchalski R.: Supporting navigators 
decisions by visualizing ship collision risk. Polish Maritime 
Research, Vol. 59, No 1, 2009, p. 83-88. 


CONTACT WITH THE AUTHOR 


Mostefa Mohamed-Seghir, Ph. D. 
Faculty of Marine Electrical Engineering, 
Gdynia Maritime University, 
Morska 81-87 
81-225 Gdynia, POLAND 
e-mail: mosem@am.gdynia.pl 


POLISH MARITIME RESEARCH, Special Issue 2012/S1 49 


POLISH MARITIME RESEARCH Special Issue S1 (74) 2012 Vol 19; pp. 50-56 
10.2478/v10012-012-0023-5 


Nonlinear observers design for multivariable 
ship motion control 
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ABSTRACT 


This paper presents the designs of two observers, which are: the extended Kalman filter 
and the nonlinear passive observer. Based on the measured values of ship position and 
heading, the observers estimate the surge, sway and yaw velocities of the ship motion. 
The observers make use of the simplified nonlinear mathematical model of ship motion in 
which the neglected ship dynamics and disturbances are modelled using bias. The designed 
observers firstly have been simulated on a computer model where their parameters were 
calibrated, and then were implemented on the physical model of the training ship “Blue 


Lady” in the ship handling centre in Ilawa-Kamionka. The comparative analysis was 
done with respect to the estimated variables describing the ship motion in three directions: surge, sway 
and yaw. 


Keywords: extended Kalman filter; nonlinear observers; ship control; dynamic positioning. 


INTRODUCTION 


Dynamic positioning is the control system which 
automatically keeps the ship position and heading using 
propellers and thrusters installed on the ship. Correctness 
of operation of the ship positioning system depends on the 
accuracy of the adopted mathematical model. Controlling ship 
motion is a difficult task, due to the fact that ship dynamics are 
nonlinear, time varying and uncertain. Moreover, environmental 
disturbances such as the wind, waves, and sea currents make the 
control task even more difficult. During a few recent decades, 
various advanced methods of ship control have been developed, 
especially in the field of robotics, which made it possible to 
reach higher quality of control. The controllers developed in this 
field require measuring all state variables, which is impractical 
and in some cases even impossible. 

In case of ship motion control in three degrees of freedom, 
the two translational velocities (longitudinal and lateral) and 
the angular speed are difficult for measuring, compared to 
the remaining states. The control systems used in dynamic 
positioning systems require the information on these velocities, 
complemented by the data on the ship position and heading. 
Usually the velocities are calculated from the position and 
heading of the ship, the measurement of which is affected 
by certain disturbances. That is why the filtration and state 
estimation play an important role in designing the dynamic 
positioning systems [1]. 

Velocity estimates are calculated from measured the ship 
position and heading data by the state observer. Unfortunately, 
these measurements are corrupted with the coloured noise 
generated by the wind, waves and oceanic currents, as well 
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as the noise generated by the sensor itself. The environmental 
disturbances acting on the ship generate two separate 
movements. The sea waves of the first order generate high- 
frequency movements, while the slowly changing forces 
generate low-frequency movements. Only slowly changing 
disturbances are to be compensated by the propellers installed 
on the ship, while the oscillatory movements generated by the 
waves (wave disturbances of the first order) should not enter 
the control system loop. This strategy is executed using wave 
filtration techniques which divide the measured position and 
heading signals into the low frequency part (LF) and the high 
frequency part generated by the waves (WF) [2]. 

The dynamic positioning systems have been developed 
since the early sixties of the last century. The first dynamic 
control systems were designed using conventional PID 
controllers working in cascade with low-pass filters or cut-off 
filters to separate the motion components connected with the 
sea waves. However, those systems introduce phase delays 
which worsen the quality of the control [1]. 

From the middle of 1970s more advanced control 
techniques started to be used, which were based on optimal 
control and the Kalman filter theory. The first solution of 
this type was presented by Balchen, Jenssen and Saelid [3]. 
It was then modified and extended by Balchen himself and 
other researchers [4, 5, 6, 7, 8, 9]. The new solutions made 
use of the linear theory, according to which the kinematic 
characteristics of the ship were to be linearized in the form of 
sets of predefined ship heading angles, with an usual resolution 
of 10 degrees. After the linearization of the nonlinear model, 
the observer based on such a model is only locally correct. This 
is the disadvantage of the Kalman filter. 


Nonlinear observers design for multivariable ship motion control 
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Fig. 1. Definition of the adopted coordinate systems: inertial frame (earth-fixed frame) xy Ze 
and body frame (body fixed frame) x,y,z, connected with the moving ship 


To avoid the linearization of the kinematic equations of ship 
motion, attempts have been made to apply nonlinear control to 
the dynamic positioning systems. Grovlen and Fossen proposed 
the use of vector backstepping as a possible solution to the 
problem of control in the dynamic positioning system. This 
solution based on the filtration of the measured ship position and 
heading signals, which were assumed to contain only the white 
noise as the disturbance. However, in practice these signals 
are also corrupted with the coloured noise connected with sea 
waves acting on the ship hull. Therefore to avoid excessive 
activity (wear and tear) of the propulsion system, the calculated 
velocity estimates were to be filtered, before introducing to 
the fedback loop, using so-called wave filtration techniques 
which extracted the wave frequency (WF) from the measured 
movements for the controller to obtain only low frequencies. 
The improved version of this design was presented by these 
same authors [10]. Those works are considered the first cases 
in which a fully nonlinear system of dynamic positioning was 
used - although this system missed wave filtration, and bias 
modelling and estimation to separate the non-modelled slowly 
changing disturbances. 

Wave filtration, bias estimation, reconstruction of LF 
movement components, and calculating noise-free estimates 
of unmeasured ship velocities were taken into account by 
Fossen and Strand who used for this purpose a nonlinear 
passive observer [11]. The nonlinear observer proposed by them 
showed new areas for controller designing, done based on the 
existing structure of the physical system and making use of the 
principle of separation developed by [12] or the backstepping 
method described by [13]. 

Activities are being in progress upon further development 
of the nonlinear systems in their application to dynamic 
positioning systems. Kim and Inman have designed a robust 
nonlinear observer which makes use of the sliding mode 
concept. The basic advantage of this observer is its robustness 
connected with the neglected nonlinearities, disturbances and 
uncertainties. The principle of operation of this observer bases 
on the work by Fossen and Strand, and it delivers the estimates 
of the linear ship velocities and the bias describing the slowly 
changing environmental disturbances [14]. Snijders, van der 
Woude, and Westhuis have designed a suboptimal observer 
making use of the Riccati equation and based on the contraction 
theory. The authors presented theoretical principles of the 
optimal nonlinear observer [15]. Hajivand and Mousavizadegan 
have designed a nonlinear observer for wave filtering and 
velocity estimation. Their observer makes use of the memory 


effect and takes into account the dependence of the added mass 
and damping coefficient changes in the ship dynamics model 
on the oscillatory forces acting on the ship hull [16]. 


MATHEMATICAL MODELLING OF THE 
SHIP 


The present work mainly aims at comparing the quality 
of operation of the Extended Kalman Filter (EKF) [17] with 
the nonlinear observer [11]. The both designs make use of the 
mathematical model of ship dynamics in which the ship motion 
is the superposition of the LF and HF movements. 


a) Kinematic equation of motion 

Two coordinate systems are introduced to describe the 
motion of a sea-going ship. The first coordinate system is 
connected with the moving ship and bears the name of the 
body-fixed frame. The centre of this system is in the centre of 
ship’s gravity. The motion of the body-fixed frame is described 
with respect to the stationary coordinate system called the 
earth-fixed frame (Fig. 1). The ship position (x, y) and the ship 
heading y, described with respect to the earth-fixed frame, can 
be expressed in the vector form as y = [x, y, y]. In the same 
form the velocities in the body-fixed frame can be given as 
v= [u, v,r]. 

The relation between the vessel-fixed and earth-fixed 
velocity vectors is given by the transformation: 


n=R@v 0) 


where R(n) is the rotation matrix defined by the formula: 


cosy -siny 0 
R(H)=R(y)=Isiny cosy 0 (2) 
0 0 1 


Notice that the matrix R(y) is non-singular for all y, hence 
R!(y) = Ry). 


b) Low frequency ship model 
For small velocities the low-frequency movements of the 
ship can be described using the following model [2]: 


Mv+Dv =1+R'()b (3) 
where t E€ R? is the vector containing the control forces and 
the moment generated by the propulsion system consisting of 
the power thrusters and water jet propellers, while the non- 
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modelled external forces and moment caused by the action of 
the wind, waves and sea currents are placed together into the 
bias component b € R’ connected with the earth-fixed frame. 

The matrix of inertia M € R”, which includes additional 
hydrodynamic inertia, is written as [2]: 


MAY 0 0 
M=| 0 m-Y, mxçş-Y; (4) 


where m is the mass of the ship, I, is the moment of inertia 
around the z axis of the coordinate system connected with the 
vessel-fixed frame, and x, is the distance between the centre 
of gravity and the origin of the coord inate system connected 
with the moving ship, X ; Y;, Y; N; and N; are the added masses 
connected with the accelerations along relevant axes. 

The linear damping matrix D € R% is defined as [2]: 


zk 2 0 
D=| 0 -Y, muy,-Y, (5) 
0 -N, MX Ug -N, 


where the velocity of motion u, = 0 in the dynamic positioning 
systems and u > 0 when the ship is in motion. Generally, 
the damping forces are nonlinear, but for the DP system the 
assumption about linear damping is satisfactory when the ship 
moves at low speed. 


c) First-order wave-induced model 

Modelling the ship movements generated by the waves of 
the first order is done based on the model of waves of the second 
order proposed by Balchen et al., [3] who used three harmonic 
oscillators. Saelid et al. improved the approximation of the WF 
wave model by introducing an additional coefficient [7]. This 
model (one for each degree of freedom) can be written as: 


i K vis 
hi,(s)=—>»* _, 6 

i s’ +20,0;8 +0; (6) 
where œ; is the dominating frequency of the wave, &; is the 
relative damping coefficient, and K , , is the parameter referring 


to the intensity of the wave. The above equation can be written 
in the state domain for three degrees of freedom as: 


Éi 03,3 13,3 . G1 4 05,3 w 
Én [+0 -2A 62 Ky oe 
Nw [0,..; bale 


Here: w, € R'is the Gaussian white noise with zero mean 
value and nų = [Kẹ Yẹ Vy’ describes the ship movements 
generated by the wave of the first order. A = diag {€,@,, 6203, 


€,@,}, Q=diag{@,, @, @,},and K=diag{K,,,, K,,, K,,4} are 
the matrices of fixed dimensions. In compact notation: 


¿= A6 +E Wy 
Nw=Cws 


(7b) 


(8a) 
(8b) 


d) Bias modelling 

It is assumed that the forces and the moment connected 
with the non-modelled dynamics and external slowly-varying 
generated by the wind, waves and sea currents are used in bias 
model. The here adopted model is the Markov process of the 
first order. 


b=-T 'b+E,w, (9) 
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where b € S? is the vector representing the slowly changing 
forces and moment, w, € R` is the vector of the Gaussian white 
noise with zero mean value, T € R™ is the diagonal matrix of 
positive time constants, and E, € R! is the diagonal matrix 
scaling the noise amplitude w,. 


e) Total system 

For conventional ships, the only available data are those 
referring to the measured ship position and heading data, 
therefore the equation relating the measured values has the 
form: 


yn tNw tY (10) 


where v € 9? is the vector of the Gaussian measurement 
noise with zero mean value. To complete these data, the ship 
observer requires measurements of u to be able to calculate 
the longitudinal and lateral forces and the torque composing 
the vector T. 

After collecting together all above equations we get the 
following model of the process: 


v=-M Dv +M'R'())b+M 1 (a) 
n=RM)v (11b) 
b=-T 'b+E,w, dlo) 
E= A6 FE, Wy (11d) 
y=n+C +v (11e) 
THE EXTENDED KALMAN FILTER (EKF) 
DESIGN 


The design of the extended Kalman filter for a dynamic 
positioning system of ship motion control will be based on the 
nonlinear model having the form: 


x = f(x) + Bu + Ew (12a) 


y=Hx+v (12b) 


where matrices f(x), B, E and H are obtained from equations 
(11a)-(11le) and take the following forms: 


R@v 03,3 
—~M'Dv +M R") M! 
f(x)= 4 B= (13a) 
=I 03,3 
A6 06,3 
03,3 033 
0 0 
E= 3x3 3x3 (13b) 
E, 033 
06.3 Ey 
H= [a 033 03,3 C] (13c) 


Here: x = [v, n, b, é]J", w = [w, w„]" and u = qt. The basic 
concept of the extended Kalman filter was proposed by Stanley 
Schmidt and sometimes this filter is referred to as the Kalman- 
Schmidt filter [18] In the obtained model the nonlinearities are 
only included in the state equation which can be written in the 
following general form: 


x(t) = f[ x(t), u(t)] + Ew(t) (14) 
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In order to obtain the practically realisable algorithm of 
state vector estimation, the nonlinear function f[x(t), u(t)] is 
expanded in the Taylor series about the current state vector 
estimate x = x: 


f[x(t), u(t)]=f[x(t), u(t)] + 


Of | x(t), u(t)] i 
+ SOO tae —x(t)] Fra 


The above expansion can be done assuming that the 
nonlinear function has partial differentials. The linear 
approximation about the current equilibrium point is obtained 
by comparing all higher-order differences to zero: 


f[x(t), u(t)]- f), u(t)] = 
_ af[ x(t), u(t)] e 
co al Bee x(t)] 


which leads to a linearized relation being the approximation 
of the equation (14): 


(15) 


(16) 


© ax = A(t)Ax + Ew(t) (17) 
where: 
Of[ x(t), u(t)] 
A(t)= — Bx | p (18) 


The extended Kalman filter is derived in the discrete form. 
For this purpose the linearized state equation (17) and the output 
equation (12b) are written in the following discrete form: 


Xp = F, Xk + GW, (19a) 


Yk = Hx, + Vv, (19b) 


where matrices F, and G, can be determined using the Euler 
forward method and the following relations: 


f 
F= Ipa +T (20a) 
OX, _ 
XKHXk 
G,=T-E (20b) 


where T > 0 is the sampling period. This system includes 
n = 15 states, the process noise covariance weight matrix 
Q= E(w! w,)e R”” and the ship position and heading 
measurement noise covariance matrix R = E (v! v,)e R`”. 
These matrices usually exist in the diagonal form. a detailed 
description of deriving the Kalman filter for the system 
described by equations (19a) and (19b) was given by Tomera 
[19]. Since in this case we deal with the exact nonlinear state 
equation (14), its discrete version will be used for calculating 
the process state values at time intervals between the sampling 
times. The measurements are done at the sampling times. 

Finally we get the following algorithm of the extended 
Kalman filter. 

The initial values: 


ar Xo (21a) 


P,- = Ef[x(0) — ¥(0)]- [x (0) — &(0)]"} = Py (210) 
Corrector: 

L, = PLH" [HPH +R] (22a) 

Ê, = ((-L,H)P, (-L,H)'+L,RL} (22b) 


x, =x, +L, (y,-Hx,) (22c) 


Predictor: 


Pk =F, PF] +G,QG, (23a) 


Xk+ =f, (Šp Up) (23b) 


where ed ee AAA is the estimated state vector, 
matrices F, and G, are determined from formulas (20a) and 
(20b), and the discrete nonlinear function F,(x,, u,) can be 
determined using the Euler forward method: 


where T > 0 is the sampling period. 


NONLINEAR OBSERVER DESIGN 


(24) 


The here presented nonlinear observer is similar to that 
proposed by Fossen and Strand and used in the dynamic 
positioning system (Fossen and Strand, 1999). The stability 
of the nonlinear observers is checked using the Lapunov 
analysis. For this purpose certain assumptions concerning the 
mathematical model of the process dynamics are to be made: 
a) The bias model (11c) and the wave model (11d) include the 
Gaussian white noises with zero mean component. In the 
Lapunov analysis these terms are to be neglected. 

Also neglected are the noises of ship position and heading 
measurements, i.e. v= 0, as these terms are of little influence 
compared to the ship movements caused by the waves. 
Itis assumed that the amplitude of the movements generated 
by the waves y, is of negligible importance and is equal 
to less than 1 degree during regular ship operation and less 
to 5 degrees in extreme weather conditions. Hence R(n) = 
=R(y)=R(y + y,). Itresults from condition (b) that R(y) ~ 
= R(y,), where y, = y + y, represents the measured ship 
heading. 

d) The following model properties are assumed for the inertia 

and damping matrices 


M=M!'>0,M=0,D>0 
After taking into account assumptions (a-d), the equations 


(11a)-(11le) describe the following mathematical model of the 
process [11]: 


b 


wm 


c 


<x 


My =-Dv+R(y)b +7 (25a) 
t =R(y,)v (25b) 
b=-T"'b (250) 
E=A,E (25d) 
y=nt+cé (25e) 


This model is used for designing the nonlinear observer. 


a) Equations of the observer 

The nonlinear observer based on the model of the ship 
and disturbances (25a)-(25e) is described by the following 
equations: 


Mv =-Dv-R'(yy)b+t+R'(yyKiy¥ (26a) 


A= R(y,)-04+ Kay (26b) 
b=-T, 'b+K3y (26c) 
E=A,E +K4y (26d) 
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j=n+C,é (260) 
where =y- ĵ € R` is the estimation error, K}, K,, KE R” 
and K} € R™ are the observer amplification matrices. 


b) Observer estimation errors 
The observer estimation errors are defined as v=v—V, 
71 =n- ñ, b=b- band é =€-€. Hence: 


MV =-DV+R'(y,)b—-R' UKY 27 
T =R) V - Kay C7) 
b =-T,'b-K3y +E,w, (270) 
E = Aye = K,y + E „Ww (270) 
y=nt+C,é (27e) 
After defining the new output: 
Z, =K,y-b=C,¥, (28) 
and vectors: 
7 
2 ~ Wo 
x, = b = | (29) 
= Wy 
E 


the error dynamics (27a)-(27e) can be written in a compact 
form as: 


M¥=—Dv+ R'(y,)C,x, = -D7 + R'(y,)z, (30a) 


Xo = A Xo +B R, N +E, w (30b) 
where: 
-K, 0; -K,C, 

A,=|-K; -T,' —-K,C, (la) 

-K, 06,3 Ay —K,C, 
C, = IK, —I,,; K,C,, | (31b) 

Lx 033 03,5 
B, =| 03,3 E, =| E, 033 (31c) 

06.3 05.3 Ey 


The observer amplification matrices can be selected in such 
a way that the error dynamics are passive. 


SIMULATION AND EXPERIMENTAL 
RESULTS 


To assess the quality and robustness of the designed 
observers, the simulation tests and experimental investigations 
were done. The experimental investigations were performed 
on the lake on which the waves are so short that they are of no 
importance for the ship motion control process, which made 
it impossible to test the quality of observers’ operation at the 
presence of waves. That is why the equations modelling the 
behaviour of waves were not included to the models of the 
observers. As a consequence, the number of the estimated state 
variables was equal to n = 9. 

The investigations were performed using the training 
ship “Blue Lady” owned by the Foundation of Sailing Safety 
and Environment Protection in Ilawa and used for training 
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navigators. The ship “Blue Lady” is an autonomous model of 
the VLCC (Very Large Crude Carrier) class tanker designed for 
transporting crude oil. The model, made of epoxide laminate 
in scale 1:24, is equipped with a set of propellers with electric 
motors fed from accumulator batteries. The overall length of 
Blue Lady is L = 13.75 m, the width is B = 2.38 m, and the 
mass is m = 22.934x 10? kg [20]. The moving coordinate system 
is fixed to the centre of gravity. 

The designed observers base on the following mathematical 
model of “Blue Lady”. 


M = diag(23664.9 24830.6 455182.1) 682) 
21.049 0 0 
D= 0 259.8 —855.4 (33) 
0 —855.4 6130.5 


The time constants for the bias were selected in the 
following way: 


T = diag(1000 1000 1000) (34) 


For the extended Kalman filter the following matrices were 
selected: 
- the process noise covariance matrix: 


Q= diag(0 0 100 10 100 100 1) (35) 


- the ship position and heading measurement noise covariance 
matrix: 


R = diag(0.01 0.01 0.01) (36) 


The amplification matrices for the nonlinear observer were 
determined as: 


K, = diag(1000 1000 1000) (37) 
K, = diag(0.5 10.5 0.025) (38) 
K, = diag(100 100 10) (39) 


Both the simulation tests and experimental investigations 
were performed using the recording frequency equal to 1 Hz. 

Firstly, the simulation tests were performed in the calculating 
environment Matlab/Simulink based on the mathematical 
model of Blue Lady described in detail by Gierusz [21, 
22]. These investigations were performed at the presence of 
measurement noises which were added to the posiotions and 
hesding measurements in simulations at the presence of the 
external disturbances. In simulation study assumed that on ship 
was acting wind with average speed equal 2 m/s in directions 
0 degrees. The simulation results are shown in Figs. 2 — 4. The 
actual and estimated velocities in surge, sway and yaw are 
shown in the bottom of plots. 

The observer algorithms which underwent simulation tests 
were then implemented on the training ship “Blue Lady” on the 
lake Silm near Ilawa. The performed experimental investigations 
aimed at assessing the quality of observers’ operation and their 
resistance to disturbances. During the investigations the ship 
position coordinates (x, y) were measured by the GPS system, 
while the ship heading y was recorded by the gyro-compass. 
The velocities: surge u, sway v and yaw r were estimated 
by the observers. Figures 5—7 show sample results of ship 
movement using three propulsion systems: the bow tunnel 
thruster (n = 0.5), the stern tunnel thruster (n = 0.5) and the 
propeller screw (n = 100 rev/min). The settings of the tunnel 
thrusters can be set within the range of <-1...+1>. 

The state variables which were estimated by the extended 
Kalman filter are shown in the left-hand columns in Figs. 6 — 7, 
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while those estimated by the nonlinear observer are collected in 
the right-hand column. The curves representing the measured 
and estimated elements of the vector yn = [x, y, y]! are smooth 
lines. The situation is different for the estimated velocity vector 
v = [u, v, r]", whose diagrams reveal certain noise in case of 
the extended Kalman filter while the estimates calculated by 
the nonlinear observer are still smooth. 
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Fig. 5. Experimental data: measured position with estimate and estimated 
velocity u in surge. Left-hand column — extended Kalman filter, right-hand 
column — nonlinear observer 
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Fig. 2. Simulation study: actual position (dotted) with estimate (solid) and 
actual velocity u in surge (dotted) and estimate (solid) 
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Fig. 3. Simulation study: actual position (dotted) with estimate (solid) and 
actual velocity v in sway (dotted) and estimate (solid) 
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Fig. 6. Experimental data: measured position with estimate and estimated 
velocity v in sway. Left-hand column — extended Kalman filter, right-hand 
column — nonlinear observer 
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Fig. 7. Experimental data: measured heading angle y with estimate and 
estimated angular rate r in yaw. Left-hand column — extended Kalman 
filter, right-hand column — nonlinear observer. When estimating the 
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Fig. 4. Simulation study: actual heading angle y (dotted) with estimate 
(solid) and actual angular rate r in yaw (dotted) with estimate (solid) 


ship motion velocity in all three degrees of freedom, we can observe the 
difference between the values estimated for the initial time, after switching 
the observers on. When starting the experimental test, the values returned 
by the nonlinear observer start from zero, while the extended Kalman filter 
detects, from the very beginning, the initial values 
of the ship motion velocity 
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REMARKS AND CONCLUSIONS 


- The article presents two algorithms of observers for 
estimating ship motion parameters in three degrees of 
freedom: surge, sway and yaw. The observers were first 
tested numerically via simulations done on a complex 
nonlinear computer model of the training ship “Blue 
Lady”. 

- The simulation tests made it possible to determine the 

parameters of the examined observers. 
Based on the satisfactory results of the simulation tests, 
the next step of investigations was the use of the examined 
observers for estimating the parameters of the real model of 
“Blue Lady” sailing on the lake Silm in Ilawa-Kamionka. 
The obtained results of the experimental investigations are 
presented in the article. Comparing the obtained results we 
can easily notice that the time-histories of the velocities 
estimated by the nonlinear observer are smooth and do not 
reveal noise, while those generated by the extended Kalman 
filter in all directions of motion reveal the measuring 
noise. 

- The designed observers can be used in dynamic multivariable 
systems controlling the motion of ships with a number of 
propeller systems, including such actions as dynamic 
positioning, tracking, or mooring. Further research 
work will be focused on designing a nonlinear dynamic 
positioning control system making use of the nonlinear 
observer. 
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Dynamic positioning system design 
for “Blue Lady”. Simulation tests 


Mirostaw Tomera, Ph. D., 
Gdynia Maritime University 


ABSTRACT 


The dynamical positioning system is a complex control consisting of a number of 
components, including: filters, observers, controllers, and propeller allocation systems. 
The design and preliminary analysis of operational quality of system operation are usually 
done based on numerical simulations performed with the aid of the mathematical model 
of the ship. The article presents a concept of the dynamic positioning system applied to 
steering the training ship Blue Lady used for training captains in the ship handling research 
and training centre owned by the Foundation for Safety of Navigation and Environment 


Protection in Ilawa/Kamionka. The simulation tests performed in the numerical environment 
of Matlab/Simulink have proved the usability of the designed system for steering a ship at low speed. 


Key words: dynamic positioning system; marine systems; ship control 


INTRODUCTION 


Automatic control of ships has been studied for over 
a century. In 1911 Elmer Sperry constructed the first automatic 
ship steering mechanism called “Metal Mike”. Today, the range 
of marine vessels covers a huge diversity of vehicles such as 
remotely operated vehicle (ROVs) and semi-submersible rigs. 
Automatic control systems for heading and depth control, way- 
point tracking control, fin and rudder-roll damping, dynamic 
positioning (DP), thruster assisted position mooring (PM) etc. 
are commercial products [1]. 

DP systems have traditionally been a low-speed application, 
where the basic DP functionality is either to keep a fixed 
position and heading or to move slowly from one location to 
another. In addition specialized tracking functions for cable 
and pipe layers, and remote operated vehicle operations have 
been available [2]. 

The dynamic positioning system for marine vessels has been 
divided into a set of dedicated modules with designed tasks. 
The most important components are shown in Fig. 1 [1]. 

The guidance system is used for planning the route for ship 
motion from the starting point to a selected target point. In the 
DP system the guidance system generates a smooth trajectory 
with position coordinates and course which lead to the next 
position. 

The unit responsible for signal processing monitors 
the measured signals and performs quality tests to detect 
extensively large variations, wild points, frozen signals, and/or 
signal drifts. Erroneous signals are to be detected and not used 
for further operations. The signal processor should check and 
evaluate the signal based on the tests of individual sensors when 


the measurements done on redundant sensors are available. 
A typical marine vessel equipped with a DP system has two 
or three gyrocompasses and the same number of position 
calculation systems. 

The main task of the observer is to provide low-frequency 
estimates of the vessel’s position, heading and velocity. Rapid 
oscillating movements caused by waves are to be filtered out. 
The observer should also be able to predict ship movements 
in the situation when ship heading and position measurements 
become unavailable (dead reckoning). 

In low-speed applications the controller calculates three 
desired parameters: the surge force, the sway force, and the 
yaw moment. Depending on the performed operation modes, 
the controller takes into account the system state estimates, the 
reference trajectory and the measured environmental conditions 
in the calculations. 

r 


Guidance Ship Thrust Actuators Environment 
System Controller Allocation ee ‘ j 
Ship dynamics 


Observer 


Signal 
Fig. 1. The major components 


of a positioning control system for marine vehicles 


LC ontrol System 


The internal logic of the controller controls the modes of 
switching between different operation types. 

The thrust allocation system maps the desired forces and 
yaw moment obtained from the controller into the required 
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propeller settings, such as the propeller rotational speed and 
pitch ratio, the angles of the rudder blade and azimuth thrusters. 
It is important that these set-points are done in an optimal 
manner, which most frequently leads to minimisation of energy 
consumption. 

The number of marine vessels with the installed DP systems 
is continuously growing due to deep-sea gas and oil mining. 
At present, the DP systems are most frequently used on shuttle 
tankers which provide services for drilling rigs. 

The beginning of the dynamic positioning systems goes 
as far to the past as to the last century’s sixties when the first 
systems acting on horizontal plane in three directions of ship 
motion: surge, sway and yaw were introduced. These systems 
made use of one-dimensional single input/single output (SISO) 
algorithms of the PID controller, along with low-pass or notch 
filters. The description of the DP systems which includes early 
stages of their development was given by Fay [3]. 

In the last century’s seventies, more advanced ship steering 
methods were introduced which based on multidimensional 
optimal control and the Kalman filter theory. The first solution 
of this type was presented by Balchen, Jenssen and Saelid [4] 
as well as by Balchen, Jenssen, Mathiasen and Saelid [5]. The 
verification of this solution on a marine ship was later presented 
by Balchen, Jenssen and Saelid [6]. Later on, this solution was 
the object of further modifications and extensions, done by 
Saelid, Jenssen and Balchen [7] who proposed a new algorithm 
of adaptation to changing frequency in order to improve the 
quality of system operation within a wide range of changes of 
environmental conditions. 

Grimble, Patton and Wise [8] presented an extended the 
analysis of the Kalman filter, comparing it to a notch filter 
which was earlier used in the dynamic positioning systems. 
Then, Fung and Grimble [9] proposed a self-tuning algorithm 
for automatic tuning of the Kalman filter matrix, obtaining good 
results. Despite the improvement in the filtration and estimation 
algorithms, the operation of the controller still based on the 
optimal control theory. 

However, some problems are observed in cases when linear 
PD controllers are used in the DP systems. Tuning the gain is 
a very complicated task which requires time-consuming tests 
done on sea with the DP system switched on. Unfortunately, 
the operational quality of the controllers changes following 
the changes of the level of environmental disturbances and 
load conditions. The DP operator has to tune manually the 
controller’s gain to adapt to changes of the environmental 
conditions. 

Another important issue is the robustness of the controller. 
The mathematical model used for modelling the ship 
motion on the sea, where the ship is subject to the action 
of wind, sea currents and waves, is strongly nonlinear and 
some phenomena are extremely difficult for mathematical 
modelling. The design of the DP controller has to take into 
account its expected robustness. Since last century’s nineties, 
this issue has been taken into account in designing linear DP 
systems by using the H,, technique. Designs of this type can 
be found in publications by Katebi, Grimble and Zhang [10], 
Kajima, Murata and Furukawa [11], Nakamura and Kajiwara 
[12], Tannuri and Donha [13], Donha and Tannuri [14] as 
well as Gierusz [15]. This type or controller meets well the 
robustness requirements in the presence of large changes of 
environmental conditions. But it is still the linear controller 
which bases on the linear model of the object; therefore 
different controllers should be designed for a number of 
operating points defined in the space of states in the vicinity 
of the point which the ship reaches during the executed 
operation. 
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In order to avoid problems connected with linearization in 
the DP systems, nonlinear controllers have been introduced. 
Fuzzy controllers were proposed by Stephens, Burnham and 
Reeve [16], Broel-Plater [17], as well as by Chang, Chen 
and Yeh [18]. Nonlinear controllers designed using the 
backstepping method and proposed by Aarset, Strand and 
Fossen [19], Strand and Fossen [20], Fossen and Grovlen [21] 
as well as by Bertin, Bittani, Meroni and Savaresi [22] were 
also successfully used. 

The publications by Fossen and Strand [23], Strand and 
Fossen [24] as well as Strand [25] present important issues of 
passive nonlinear observers with adaptive wave filtration. An 
advantage of the use of the nonlinear theory of passiveness was 
the reduced complexity of the code packages used for steering. 
Pettersen and Fossen [26], Pettersen, Mazenc and Nijmeijer 
[27] as well as Bertin, Bittani, Meroni, and Savaresi [22] have 
worked out the DP control algorithms for seagoing vessels in 
which the number of propellers is smaller than the number of 
freedom degrees (so-called under-actuated vessels). Agostinho, 
Moratelli, Tannuri and Morishita [28] as well as Tannuri, 
Agostinho, Morishita and Moratelli [29] have proposed the 
use of nonlinear sliding control in the DP system. 

The application of hybrid control theory proposed by 
Hespanha [30], Hespanha and Morse [31], Hespanha, Morse 
and Liberzon [32], as well as the fault-tolerant control proposed 
by Blanke, Kinnaert, Lunze and Staroswiecki [33] has made it 
possible to design a correct control architecture for integrating 
multifunction controllers which link discrete events and 
continuous steering. The effect of operation of such a controller 
applied in DP systems was described by Sorensen, Quek and 
Nguen [34], Nguyen [35], Nguyen and Sorensen [36] who 
proposed a design of the controller with the superior switching 
logic to select a proper controller and observer from a set of 
controller and observers assumed for different environmental 
conditions. 

The applicability of the DP systems on shuttle ships, where 
such operations as position keeping, sailing along a given 
trajectory, and unloading with FPSO (floating production 
storage and offloading) are performed, was analysed by 
Morishita and Cornet [37], Morishita, Tannuri and Bravin 
[38], Tannuri and Morishita [39], as well as Tannuri, Saad and 
Morishita [40]. 

The article presents a design of the DP system applied 
to steering the training ship Blue Lady. The investigations 
performed on the observers which could be hypothetically 
used in the designed DP system had been done earlier and 
were presented in the following articles and papers: the discrete 
Kalman filter [41], the continuous Kalman-Bucy filter [42], and 
the extended Kalman filter and nonlinear observer [43]. 


THE MATHEMATICAL MODEL OF BLUE 
LADY DYNAMICS FOR LOW SPEED 


The training ship Blue Lady is owned by the Foundation for 
Safety of Navigation and Environment Protection in Ilawa and 
is used for training captains to perform complex and difficult 
manoeuvres of a large ship. The ship, made of epoxide laminate 
in the scale of 1:70, is a replica of a tanker used for transporting 
crude oil. The overall length and breadth of the physical model 
are, respectively, Loa = 13.75 [m], and B = 2.38 m. In the full- 
load state its mass is m = 22830 [kg], moment of inertia I, = 
436830.3 [kgm?], while center of gravity is located at midship, 
hence x, = 0 [m]. The model is equipped with a set of actuators 
with electric motors fed from an accumulator battery, including 
the main propeller with a conventional plane rudder, and four 
jet thrusters: two tunnel (the bow thruster and the stern thruster) 
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and two rotational jet propellers with changing rotational speed 
(the bow propeller and the stern propeller). The distribution of 
actuators on Blue Lady is presented in Fig. 2. 

The presented article takes into account only three actuators, 
namely the main propeller and two tunnel thrusters: bow and 
stern. 

A complex mathematical model of Blue Lady dynamics, 
complemented by the modelled actuators which relatively well 
model its real behaviour, was developed by Gierusz [44]. 


a) Simplified mathematical model 


The mathematical model of ship dynamics is described by 
the position vector ņ = [x, y, y] which comprises the position 
coordinates: North x, East y, and the heading y. These quantities 
are calculated in the Earth-fixed coordinates fixed to the water 
region map (n-frame). The velocity components: surge u 
and sway v as well as the yaw rate r form the velocity vector 
v= [x, v, r] calculated into the body frame fixed to the moving 
ship (b-frame) and hence moving along with it. 


Fig. 2. Distribution of actuators on Blue Lady 


The dynamic positioning system is designed for steering the 
ship at low speed. The proposed multidimensional controller 
applied in the DP system is determined using the mathematical 
model of the ship sailing at low speed, given by the following 
formulas [45]: 


=Rv (1) 


Mv+D,v =1 (2) 
where R(y) is the matrix of rotation calculated from the 
formula: 


cosy -siny 0 
R(y) =| siny cosy 0 (3) 
0 0 1 


The mass matrix M comprises the parameters of inertia of 

a rigid body, its dimensions, weight, mass distribution, volume, 
etc., and the added mass coefficients: 

m-X, 0 0 
M= 0 m=Y, mk.=-Y: (4) 
0 mrta-N; L—N. 
The linear damping matrix D, is connected with the 
hydrodynamic damping forces and is calculated for a selected 


constant small value of the surge speed v = vy = [Up, 0, 0] 
[46]. 


=x 0 0 
Di=| 0 -Y, =] (5) 
O -N =. 


It is assumed that the both vectors ņ and v are measured. 
The parameters of the Blue Lady model motion calculated for 
the velocity uy = 0.1 [m/s] are given in Tab. 1. 


Tab. 1. Calculated values of parameters of the simplified mathematical 
model of Blue Lady 


Parameter 


Parameter 


730.5 
1896.2 
18351.9 
0 


b) Simplified models of propellers 


This section describes the mathematical models of 
propellers which were used in ship steering. Their operation 
is analysed for low speed ranges. The thrust force for the 
propeller screw is: 


F\=k Joo, (6) 
where: 
k; = 1.9589. 


The thrust forces for the tunnel thrusters: bow F, and stern 
F, are given by the formulas: 


F;=k@,ie {2,3} (7) 


where: 
k, =k, = 44145. 


The surge and sway forces generated by those thrusters: 
u = [F FF)" (8) 


The vector of the forces acting on the ship in relation to the 
rotational speed of the propellers œ;: 


t= Tu (9) 
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Tx 1 0 0 F (œ) 
ty |=|0 1 1 F,(@,) (10) 
TN 0 L, —-L;] LK@;) 

where: 

L, = 3.24 [m], 

L, = 2.376 [m]. 


The matrix T is the thrusters configuration matrix. 
GUIDANCE SYSTEM 


The ship guidance system was designed for performing 

two tasks: 

- Point stabilisation. The task of the guidance is to keep the 
vessel at the given point and constant heading; 

- Trajectory tracking. The task consists in moving the ship 
along the time-parameterised reference trajectory. 


The mathematical model of the ship moving in the 
horizontal plane describes the surge dynamics, the sway 
dynamics and the yaw dynamics and is the model with three 
degrees of freedom (3 DOF). The description of the ship 
motion in the dynamic positioning system in which the ship 
keeps a fixed position or moves slowly from one position to 
another is done in a very small area and there are no problems 
with mapping the spherical shape of the Earth onto the plane. 
Fig. 3 shows schematically an example motion of the vessel 
from the starting point P to the final point K. This motion 
was defined using three coordinate systems. The first system 
is the Cartesian coordinate system x"y" fixed to the map of 
the water region and represents the plane tangential to the 
surface of the Earth in the region in which the manoeuvre is 
performed. In this system the x" axis is directed north, while 
the y” axis is directed east. This reference system is related 
to the navigation on the Earth surface and for simplicity 


(North) 


x” 


Ya 


bears the name of the n-frame. The second coordinate system 
xby> is the relative system fixed to the moving ship. In this 
system the x? axis is directed towards the longitudinal axis of 
the ship (from stern to bow), while the y? axis is the lateral 
axis directed towards the starboard side. The ship position 
(x, y) and heading y are calculated with respect to the absolute 
coordinate system (n-frame) and collected in the position 
vector q = [x, y, w]', while the linear velocity components 
(u, v) and the angular speed r are calculated in the relative 
coordinate system xy and collected in the velocity vector 
= [x, v, r]'. The origin of the relative coordinate system is 
usually situated in the centre of ship gravity. Moreover, the 
description of the motion of the ship changing the position 
makes use of the third coordinate system x'y' playing the role 
of the reference coordinate system and bearing the name of 
the r-frame. This system is used for describing the ship motion 
from the starting point P to the final point K. It has the same 
properties as the n-frame, the only difference is that the x" 
axis changes direction and the coordinates of the origin of 
this system are shifted to the starting point P of the currently 
passed trajectory segment. 

Fig. 3 shows an example situation in which the ship 
situated at the starting point P has the position coordinates 
np = p> Yp» welt and has to change the position moving to 
the final point K, where it will have the position coordinates 
N= [Xe Yeo Yl. These position coordinates are written in 
the n-frame. 

During position changes, the movements executed by 
the ship are governed by the DP controller, the operation of 
which consists in such ship steering that the ship follows with 
a satisfying accuracy the positions of the moving reference 
system x‘y4, fixed to the virtual ship. All this motion from 
the starting point P to the final point K is described in the 
r-frame. 

To do this, the ship position coordinates at the final point 
K n= Rk Ye Wil, the coordinates of the current ship 
position ņ = [x, y, w]’ and the required ship position 


n 


y 
n (East) 


Fig. 3. Definitions of the introduced coordinate systems and locations of the vessel changing the position 
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1,=[Xy. yi Wy]! are recalculated to the r-frame. The starting 
coordinates and the orientation of the r-frame are identical with 
the ship position at the starting point P. 

The coordinates of the final point toward which the ship 
moves are calculated in the r-frame from the following 
relation: 


nk =k -nP): RO) (11) 
where: 
R(@) -— the rotation matrix defined in the following way: 
cosd —singd 0 
R@)=| sind cosd 0 (12) 
0 0 1 
while the angle of rotation: 
b=-Wp (13) 


Let ‘rf =R}, yi J" be the required position of the ship and 
yi, the required course written in the r-frame. The required 
ship position and course in the r-frame, which control the 
ship movement from point P to point K is given by the vector 
nL (= r t) y] , shortly written: 


Na = Nea (14) 


The components of the position 7,,(t), velocity 1,,(t) and 
acceleration N4 (t) of the virtual ship moving from the starting 
point P to the final point K are calculated in the reference frame 
shown in Fig. 4. The required input values of the velocity 1,,(t) 
and the acceleration 1,,(t) along the trajectory should not exceed 
physical limits valid for the ship; therefore relevant limits for 
the signals were introduced to the system shown in Fig. 4. 

The desired velocities calculated in the system presented in 
Fig. 4. required recalculation from fixed reference coordinate 
system r-frame to the body frame fixed to the moving vessel 
(b-frame) according the following formulas: 


ta =v = R yenii 
fa =Va = VSR OA -R' We) Ña (16) 


where: 
Ye = Y — Ya the matrix of rotation R(y,) is calculated from 
equation (3), while the matrix S has form 


(15) 


0 -1 0 
S=|1 0 0 (17) 
0 0 0 
tra) 


ya (t) 


Tira È) 
> 
Fig. 4. Reference frame for generating input signals for the DP controller 


FULL STATE FEEDBACK PD TRACKING 
CONTROLLER 


The goal of the control is to track the smooth trajectory 
(n,(t).n,(t).1,(t)) generated by the ship guidance system. 
The control system calculates the deviations from the set input 
values. The position error n, is calculated in the absolute r- 
frame, while the velocity h, and the acceleration Ñ, errors are 
calculated in the relative b-frame fixed to the moving ship. 


Ne =N—Na (18) 
Ne =n -Na (19) 
Ne =n -Na (20) 


The DP control system includes a controller which stabilises 
the error dynamics, i.e. 


H+ Kn. + Kpn. =0 (21) 


After converting the equation (21) to get the ship 
acceleration 1 we arrive at: 
N= Na-Kpn.- Kn. (22) 
Based on the equation (2), the proposed rule of control is 
described by the formula: 
t = My + Dyv (23) 


The velocity vector derivative v in the equation (23) is 
determined from the kinetic equations (1) differentiated with 
respect to time, which leads to: 


v=R' Qi- R] 24) 


The derivative of the rotation matrix R(y) is determined 
from the formula: 


R(y) = rR(y)S (25) 
where: 

r=W. 

After placing the relation (22) into equation (24), and then 
into the proposed rule of control (23), we arrive at the formula 
describing the algorithm of operation of the proposed DP 
controller: 


t= MR (tig - Kpr,- Kpn- R(W)v] + Diy 26) 
CONTROL ALLOCATION 
The problem with control allocation can appear when the 


number of actuators is larger than the number of the controlled 
degrees of freedom. The algorithm can be realised in two steps, 


as shown in Fig. 5. 
Inverse @, 
Mapping 


Fig. 5. The thrust allocation problem 


Te Force 


Allocation 


In the first step, bearing the name of force allocation, the 
desired generalised force t, is decomposed into all propellers 
under consideration. The quality of the decisions made in this 
step depends on how good the force allocation algorithm is. 
The second step consists in finding such settings in the actuators 
which generate the required forces F. This step is referred to as 
the inverse transformation as it consists in finding the inverse 
characteristics of the devices. 

In the examined DP system three propellers are used for 
controlling the ship motion. Their configuration is given in 
Fig. 2. In this case the force allocation consists in solving the 
equation: 


u=T't (27) 
E] fl 0 oO] [ty 
FRl={0 1 1 | fry (28) 
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Since the number of actuators is equal to the number of 
the controlled degrees of freedom, there is no problem with 
calculating the matrix T-!. In case the number of actuators is 
greater than the number of the number of the controller degrees 
of freedom, the propeller configuration matrix T~! does not 
exist and we cannot find a direct solution. In this case the 
optimisation methods are to be used [47]. 

After decomposing the forces into particular propellers, 
from the transformed formulas (6) and (7) we obtain the desired 
values for the main screw propeller: 


o, =sen(F) )/R /k, (29) 
and for the tunnel thrusters: 
@=F,/k;, i¢{2,3} (30) 


The desired rotational speed for the main propeller is limited 
to +200 [rev/min], while for the bow and stern tunnel thrusters 
the limits are within the range of +1 [-]. 


NUMERICAL SIMULATION RESULTS 


The simulation tests were performed in the mathematical 
environment Matlab/Simulink. Particular components of the 
control system shown in Fig. 1 were modelled as blocks, 
while the algorithms describing the action of those blocks 
have been written in the form of S-functions in the Matlab 
code. Finally, the algorithms of the DP control system block 
components were translated into S-functions written in the 
code C++, 

The operational quality of the designed positioning control 
system was tested on a complex mathematical model of the 
training ship Blue Lady worked out by Gierusz [15, 44]. Based 
on the complex mathematical model, a simplified model of ship 
dynamics was worked out and then used for the synthesis of the 
controller (26). In the simulation tests the following parameters 
of the controller were assumed: 


K, = diag(1 1 1) (31) 
K, = diag(100 100 100) (32) 


The values of the parameters in the superior system were 
the following: 


Z = diag(1 1 1) 
Q= diag(0.012 0.01 0.015) (34) 


The DP controller applied for determining the desired forces, 
qT, requires the information about six state variables describing 
the motion of the ship, out of which the coordinates of ship 
position (x, y) and course y are measured and collected in the 
position vector n = [x, y, y], while the velocity components: 
surge u, and sway v, as well as the yaw rate r, which are not 
measured, are estimated by the nonlinear observer described 
in detail by Tomera [43] and collected in the velocity vector 
v=[kx, v, r]. 

An example case of ship motion at low speed, shown in 
Fig. 6, consisted in moving the ship from one quay to another. 
In that time the ship got away from the first quay moving in 
the lateral direction and stopped at the first stopping point. 
Then it sailed in longitudinal direction to the other stopping 
point at which it performed the rotating manoeuvre. When 
it took a proper position with respect to the second quay, 
the ship sailed astern along a distance and finally performed 
the manoeuvre of touching the land at the second quay. The 
entire trajectory covered by the ship was divided into five 
manoeuvres. Each manoeuvre consisted in changing the 


(33) 
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parameters of the ship position vector ņ from one steady 
state to another. 

The set values for the changing components of the ship 
position vector n were generated in the guidance system. The 
recorded time-histories of particular vector components are 
shown in Fig. 7. The ship position coordinates n are expressed 
in the fixed coordinate system, the r-frame, moved to a new 
position for each manoeuvre. 
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Fig. 6. Simulation results - example ship motion trajectory 
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Fig. 7. Simulation results — the solid lines represent positions while the 
dashed lines represent reference (r-frame) n, = x - position x, 
N, =y - position y, n; = y - heading 


In transient states, when the ship moved from one fixed 
point to another, its motion was executed at some speed which 
was also set by the guidance system. The recorded time- 
histories of the velocity components: both the set ones and those 
executed by the ship are shown in Fig. 8. These velocities are 
expressed in the relative coordinate system, the b-frame, fixed 
to the moving ship. 

The components of the desired forces which are determined 
by the DP controller and then should be executed by the 
propellers are shown in Fig. 10. The ship’s DP controller 
calculated these forces based on the errors between the desired 
and measured values of the position and velocity vectors and 
the desired acceleration vector values (Fig. 8). 

In the force allocation system the desired values of the 
force vector t were recalculated to the desired parameters of 
the thrusters used for the steering. The time-histories of the 
desired values calculated for the thrusters are shown in Fig. 11. 
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Fig. 8. Simulation results — the solid lines represent velocities while the 
dashed lines represent references (b-frame) v, = u - surge, v, = V - sway, 
v3 5r - yaw rate 
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Fig. 9. Simulation results — the doted lines represent set accelerations 
(b-frame) v,= ù - surge, ¥2=¥ - sway, V3 =F - yaw 
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Fig. 10. Simulation results — virtual forces generated by the control law 
T, - surge force, T, - sway force, T; - yaw moment 


The figure does not include current values of these thrusters, 
because they cannot be measured in the real system. The 
rotational velocity œ| of the main screw propeller is expressed 
in revolutions per minute, while the desired rotational velocities 
©, and œ; for the bow and stern tunnel thrusters, respectively, 
are expressed in the normalised numbers ranging between 
[1, +1]. 
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Fig. 11. Simulation results — computed propeller angular velocities by the 
allocation algorithm œ, - desired rotational velocity of the main propeller, 
@, - desired rotational velocity of the tunnel thruster mounted at the bow, 
w; - desired rotational velocity of the tunnel thruster mounted at the stern 
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CONCLUSIONS 


The article presents a general concept and the design of the 
ship motion control at low speed. The operational quality of the 
designed system was tested using simulation calculations. The 
training ship “Blue Lady” was selected as the object of control, 
as it provides opportunities for future experiments in the ship 
handling research and training centre in Ilawa Kamionka 
to verify the results of the simulation tests. The performed 
simulation tests did not take into account the environmental 
disturbances which normally act on the ship. The planned future 
experiments will make it possible to determine the level of 
disturbances acting on the training ship “Blue Lady” when it 
sails on the lake and take them into account in the mathematical 
model of the control system. 

The article is a first proposal of the multidimensional 
control algorithm applied to steering the training ship ,,Blue 
Lady”, which in the future will be improved by introducing the 
effect of environmental disturbances. Other analysed directions 
of system development include the use of the multidimensional 
PID controller instead of the presently used PD controller 
and the development of the propeller allocation system by 
including the rotational jet propellers installed on the training 
ship „Blue Lady”. The presently used propeller allocation 
system does not require optimisation as the number of the 
used propeller is equal to the number of degrees of freedom 
in which the steering takes place. The use of two additional 
thrusters and jet propellers: one at the bow and one at the stern, 
will introduce four additional set values, as each additional 
rotational jet propeller requires the set values of the setting 
angle and the rotational speed. 


BIBLIOGRAPHY 


1. Lindegaard K.-P.: Acceleration Feedback in Dynamic 
Positioning, PhD Thesis, Norwegian University Science & 
Technology, Deptment of Engineering Cybernetics, Trondheim, 
Norway, 2003 

2. Sorensen A.J.: a survey of dynamic positioning control systems. 
Annual Reviews in Control, Vol. 35, Issue 1, pp. 123-136, 
Elsevier Ltd., 2011 

3. Fay H.: Dynamic positioning systems, principles, design and 
applications. Editions Technip, Paris, France, 1990 

4. Balchen J.G., Jenssen N.A., Saelid S.: Dynamic positioning 
using Kalman filtering and optimal control theory. In JFAC/IFIP 


POLISH MARITIME RESEARCH, Special Issue 2012/S1 63 


Dynamic positioning system design for “Blue Lady”. Simulation tests 


Symposium on Automation in Offshore Oil Field Operation. 
Amsterdam, The Netherlands, pp. 183-186, 1976 

5. Balchen J.G., Jenssen N.A., Mathiasen E., Saelid S.: a dynamic 
positioning system based on Kalman filtering and optimal 
control, Modeling, Identification and Control, Vol. 1, No 3, pp. 
135-163, 1980 

6. Balchen J.G., Jenssen N.A., Saelid S.: Dynamic positioning of 
floating vessels based on Kalman filtering and optimal control, 
Proceedings of the 19* IEEE conference on Decision and 
Control, New York, pp. 852-864, 1980 

7. Saelid S., Jenssen N.A., Balchen J.G.: Design and analysis of 
a dynamic positioning system based on Kalman filtering and 
optimal control, EEE Transactions on Automatic Control, Vol. 
28, No 3, pp. 331-339, 1983 

8. Grimble M.J., Patton R.J., Wise D.A.: Use of Kalman Filtering 
Techniques in Dynamic Positioning Systems. JEE Proceedings, 
Vol. 127, Part D, No 3, pp. 93-102, 1980 

9. Fung P. T-K., Grimble M.: Dynamic ship positioning using self- 
tuning Kalman filter. JEEE Transactions on Automatic Control, 
Vol. 28, Issue 3, pp. 339-349, 1983 

10.Katebi M.R., Grimble M.J., Zhang Y.: H,, robust control design 
for dynamic ship positioning. JEE Proceedings - Control Theory 
and Applications, Vol. 144, No. 2, pp. 110-120, 1997 

11.Kijima K., Murata K., Furukawa Y.: Design of the dynamic 
positioning system for self-propulsive barge under external 
disturbances, Proceedings of the IFAC conference on control 
applications in marine systems, (CAMS 98), Fukuoka, Japan., 
pp. 69-76, 1998 

12.Nakamura M., Kajiwara H.: Control system design and model 
experiments on thruster assisted mooring system. Proceedings of 
the 7" international offshore and polar engineering conference 
(ISOPE), pp. 641-648, 1997 

13.Tannuri E.A., Donha D.C.: H,, controller design for dynamic 
positioning of turred moored FPSO. Proceedings of the 5" 
IFAC conference on manoeuvering and control of marine crafts, 
(MCMC-2000), Aalborg, Denmark. pp. 269-275, 2000 

14.Dohna D.C., Tannuri E.A.: Non-linear semi-submersible 
positioning system design using an H,, controller. Proceedings 
of the 5 IFAC Conference on Control Application in Marine 
Systems(CAMS-2001), Glasgow, Scotland, 2001 

15.Gierusz W., Synthesis of multivariable systems of precise ship 
motion control with the aid of selected resistant system design 
methods, Gdynia Maritime University Publishers, (in Polish), 
2005 

16.Stephens R.I., Burnham K.J., Reeve P.J.: a practical approach 
to the design of fuzzy controllers with application to dynamic 
ship positioning. Proceedings of the IFAC Conference on control 
applications and marine systems, (CAMS-95), Trondheim, 
Norway, 1995 

17.Broel-Plater B.: Fuzzy control for vessel motion, Proceedings 
of the 4" International Symposium on Method and Models in 
Automation and Robotics, Vol. 2, pp. 697-702, Miedzyzdroje, 
Poland, 1998 

18.Chang W.-J., Chen G.-J., Yeh Y.-L.: Fuzzy control of dynamic 
positioning systems for ships, Journal of Marine Science and 
Technology, Vol. 10, No. 1, pp. 47-53, 2002 

19.Aarset M.F., Strand J.P. Fossen T.I.: Nonlinear vectorial observer 
backstepping with integral action and wave filtering for ships, 
Proceedings of IFAC Conference on Control Applications in 
Marine Systems, (CAMS-98), Fukuoka, Japan, pp. 83-88, 1998 

20.Strand J.P., Fossen T.I.: Nonlinear Output Feedback and 
Locally Optimal Control of Dynamically Possitioned Ships: 
Experimental Results, In: Procedings of IFAC Conference on 
Control Applications in Marine Systems (CAMS-98), Fukuoka, 
Japan, pp. 89-94, 1998 

21.Fossen T.I., Grovlen A.: Nonlinear output feedback control of 
dynamic positioned ships using vectorial observer backstepping, 
IEEE Transactions on Control Systems Technology, Vol. 6, No. 
1, pp. 121-128, 1998 

22.Bertin D., Bittani S., Meroni S., Savaresi S.M.: Dynamic 
positioning of a single-thruster vessel by feedback linearization, 
Proceedings of the 5* IFAC conference on manoeuvering and 


64 POLISH MARITIME RESEARCH, Special Issue 2012/S1 


control of marine crafts, (MCMC-2000), Aalborg, Denmark. 
pp.281-286, 2000 

23.Fossen T.I., Strand J.P.: Passive nonlinear observer design 
for ships using Lyapunov methods: Experimental results with 
a supply vessel, Automatica, Vol. 35, No. 1, pp. 3-16, 1999 

24.Strand J.P., Fossen T.I.: Nonlinear passive observer design for 
ships with adaptive wave filtering. New Directions in Nonlinear 
Observer Design (Nijmejer H., Fossen T.I., Editors). pp. 113- 
134. Springer-Verlag, London Ltd., 1999 

25.Strand J.P.: Nonlinear position control systems design for 
marine vessels, PhD Thesis, Norwegian University Science & 
Technology, Deptment of Engineering Cybernetics, Trondheim, 
Norway, 1999 

26.Pettersen K.Y., Fossen T.I.: Underactuated Dynamic Positioning 
of a Ship - Experimental Results, JEEE Transactions on Control 
Systems Technology, Vol. 8, No. 5, pp. 856-863, 2000 

27.Pettersen K.Y., Mazenc F., Nijmeijer H.: Global uniform 
Asymptotic Stabilization of an Underactuated Surface Vessel: 
Experimental Results, Transactions on Control Systems 
Technology, Vol. 12, No. 6, pp. 891-903, 2004 

28.Agostinho A.C., Moratelli L., Tannuri E.A., Morishita H.M.: 
Sliding mode control applied to offshore dynamic positioning 
systems, Proceedings of the 8 IFAC International Conference 
on Maneouvering and Control of Marine Craft, MCMC-2009, 
Guaruja Brazil, 2009 

29.Tannuri E.A., Agostinho A.C., Morishita H.M., Moratelli L.: 
Dynamic positioning systems: An experimental analysis of 
sliding mode control, Control Enginering Practice, Vol. 18, No. 
10, pp. 1121-1132, 2010 

30. Hespanha J.P.: Tutorial on Supervisory Control, Lecture 
Notes for the Worhshop Control using Logic Switching for the 
40" Conference on Decision and Control, Orlando Florida, 
2001. 

31.Hespanha J.P., Morse A.S.: Switching between Stabilizing 
Controllers. Automatica, Vol. 38, No. 11, pp. 1905-1917, 2002 

32.Hespanha J.P., Liberzon D., Morse A.S.: Hysteresis-based 
Switching Algorithms for Supervisory Control of Uncertain 
Systems. Automatica, Vol. 39, No 2, pp. 263-272, 2003 

33.Blanke M., Kinnaert M., Lunze J., Staroswiecki M.: Diagnostics 
and Fault-Tolerant Control. Springer-Verlag, Berlin, Germany, 
2003 

34.Sorensen A.J., Queck S.T., Nguyen T.D.: Improved operability 
and safety of DP vesels using hybrid control concept. 
Proceedings of the International conference on technology & 
operation of offshore support vessels (OSV), 20 - 21 September, 
Singapore, 2005 

35.Nguyen T.D.: Design of Hybrid Marine Control Systems 
Jor Dynamic Positioning, PhD Thesis, Department of Civil 
Engineering, National University of Singapore, 2006 

36.Nguyen T.D., Sorensen A.J.: Switching control for thruster- 
assisted position mooring. Control Enginering Practice, Vol. 17, 
No. 9, pp. 985-994, 2009 

37.Morishita H.M., Cornet B.J.J.: Dynamics of a turret-FPSO and 
shuttle vessel due to current. Proceedings of IFAC Conference 
on Control Applications in Marine Systems, (CAMS-98), 27-30 
October, Fukuoka, Japan, pp. 93-98., 1998 

38.Morishita H.M., Tannuri E.A., Bravin T.T.: Methodology for 
dynamic analysis of offloading operations, Proceedings of 
IFAC Conference on Control Applications in Marine Systems, 
(CAMS-2004), July 7-9, Ancona, Italy, pp. 459-464, 2004 

39.Tannuri E.A., Morishita H.M.: Experimental and numerical 
evaluation of a typical dynamic positioning system, Applied 
Ocean Research, Vol. 28, No. 2, pp. 133-146, 2006 

40.Tannuri E.A., Saad A.C., Morishita H.M.: Offloading operation 
with a DP shutle tanker: Comparison between full scale 
measurement and numerical simulation results. Proceedings of 
the IFAC International Conference on Maneuvering, Guaruja, 
Brazil, 2009 

41.Tomera M.: Discrete Kalman filter design for multivariable 
ship motion control: experimental results with training ship, 
Joint Proceedings of Gdynia Maritime Academy & Hochschule 
Bremerhaven, Bremerhaven, pp. 26-34, 2010 


Dynamic positioning system design for “Blue Lady”. Simulation tests 


42.Tomera M.: Kalman-Bucy filter design for multivariable ship 
motion control, Methods and Algorithms in Navigation — Marine 
Navigation and Safety of Sea Transportation, Editors Adam 
Weintrit & Tomasz Neumann, Published by CRC Press/Balkema, 
pp. 21-31, 2011 

43.Tomera M., Nonlinear observers design for multivariable ship 
motion control, Polish Maritime Research, Accepted for print. 

44.Gierusz W.: Simulation model of the shiphandling training 
boat „Blue Lady”, Proceedings of the 5" IFAC Conference on 
Control Application in Marine Systems (CAMS-2001), 18-20 
July, Glasgow, Scotland, 2001 

45.Fossen T.I.: Guidance and control of ocean vehicles, John Wiley 
& Sons, Chichester, UK, 1994 

46.Fossen T.I.: Marine Control Systems: Guidance, Navigation, 
and Control of Ships, Rigs and Underwater Vehicles, Marine 
Cybernetics, Trondheim, Norway, 2002 


47.Fossen T.I.: a Survey of Control Allocation Methods for 
Ships and Underwater Vehicles, Proceedings of the 14" IEEE 


Mediterriaen Conference on Control and Automation, June 28- 


30, Ancona, Italy, 2006 


CONTACT WITH THE AUTHOR 


Mirostaw Tomera, Ph. D. 
Faculty of Marine Electrical Engineering, 
Gdynia Maritime University, 
Morska 81-87 
81-225 Gdynia, POLAND 
e-mail: tomera@am.gdynia.pl 


POLISH MARITIME RESEARCH, Special Issue 2012/S1 


65 


POLISH MARITIME RESEARCH Special Issue S1 (74) 2012 Vol 19; pp. 66-72 
10.2478/v10012-012-0025-3 


Measuring system for parallel moving ships 


Anna Waszkiel, M.Sc., 
Gdynia Maritime University 


ABSTRACT 


The paper introduces algorithm for determining the relative positions of two ships 
manoeuvring as a pair. This algorithm also takes into account determination of angle y, 
which is difference between present approaching vessel and guidance vessel course. Relative 
positioning system is a vision system based on three colours LEDs matrix and rotating 
CCD camera. There are presented ways of distance calculation based on photogrammetric 
methods from the known distance between the characteristic points of the real. Several 
possible cases are taken into account. The considerations have been illustrated on the 


basis of model of the system. These results confirm the correctness of the operation of the algorithm that 
is used by the designed measuring system. 


Keywords: vision system; relative ship position determination; safe ship control; computer simulation 


INTRODUCTION 


Each ship manoeuvring on the marine area, for safety 
navigation, should have information about courses and 
position of the other vessels navigating in the vicinity. Safety 
of units being in motion is subjected to the consideration of 
the hydro-meteorological factors, navigational dangers (eg. 
shallow waters, rocks, wrecks) and other vessels positions in 
multidimensional steering process. In ships control system 
information form electro navigational systems like GPS, AIS, 
radar and ship’s log is used for this purpose. This information 
can be derived also by navigator’s decision support system, 
like ARPA. Taking into account two vessels moving in parallel, 
safety of the manoeuver depends not only on proper ships 
location in navigational space, but also on accurate relative 
position determination. Underway Replenishment Operations 
(UNREP), fishing boats trawling in pair, pipe laying operations 
and some research situations are real examples of parallel ships 
movement on the sea. 

For the safe performance of manoeuver in vessels parallel 
motion, it is necessary to pinpoint transversal and longitudinal 
distance between ships and also to determine parallelism of 
the course. In case of no parallelism of ship’s courses it is 
necessary to determine angle y, which is difference between 
present approaching vessel and guide ships course, which is 
shown in the Fig. 1. 

When taking into account manoeuver of ships parallel 
movement, whole procedure consisting of three phases should 
be taken into consideration. It is not only case when vessels 
are arranged side to side amidships. This exemplary operation 
is divided into three phases: approaching, parallel movement 
and departure, which is shown in the Fig. 2. 
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Fig. 1. Approaching and guide ships relative position coordinates 
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Approaching and departure are the most dangerous phases, 
because of big influence of forces and yawing moments acting 
on vessels bodies. Parallel movement is easier to control, 
because only suction force is significant, there is hardly any 
yaw moment there due to suction phenomenon [7]. So it is 
important to develop relative position counting algorithm, 
which will work properly not only in parallel movement, but 
also in first and third phase. 


Measuring system for parallel moving ships 


METHODS FOR DETERMINING 
RELATIVE POSITIONS OF THE SHIPS 


In case of two ships manoeuvring on the sea surface only 
contactless measurement methods can be applied. Distance 
can be measured directly or computed analytically from shape 
measurement on the basis of mathematical transformations. 
This is possible using sound waves (ultrasonic rangefinders) 
[10] or light waves [2, 5] (optical methods) as a medium. 
Relative position and course determination in marine navigation 
systems is based on radar data. Radar angular discrimination 
should be less than 1° and distance discrimination should reach 
20m or 1.5% of range according to IMO Performance Standards 
[6]. Taking radar into account there is also problem of shadow 
areas and reflections from water surface in the vicinity of the 
ship. Discrimination values for ships models are to big to rely 
on them in relative position estimation system. In addition 
radar in its working principle uses microwave radiation, which 
is unhealthy for the crew working in a small distance from the 
other ship. Radar could be treated as additional element of the 
whole positioning system during approaching and departure 
phase. 

Optical measurement methods are based on measuring the 
reflection or scattering on the surface of the test object. The 
most common method of the optical measurement is based on 
the laser triangulation. Triangulation phenomenon was first 
time mathematicaly described by W. Snell van Royen in 1615. 
Over the past 20 years there has been seen development of this 
method [5], caused by CCD cameras widespread usage and 
laser technology development. The next well known optical 
distance measurement method is structured light projection. The 
method is based on raster patterns projection onto the surface of 
an object. Then observations carried out with a camera system 
are analysed. Received patterns are processed and on the basis 
of analytical dependences the distance is obtained [4]. 

For a mutual vessel’s position evaluation system it is not 
possible to apply any of the presented above methods. Cone- 
shaped bow makes laser range-finders distance measurements 
unusable. There is no simple way to verify to the which side 
of the vessel the measurement was made, when the y angle is 
negative and/or big (Fig. 1). Furthermore computer analysis 
would require knowing the exact shape of the hull. This would 
make the system possible to use only on board of certain type 
of ships. Structural light projection is a method which for 
physical reasons — like strong solar light, dark scratched hull, 
could not give desired effects. Of these reasons decision was 
made, to use photogrammetric methods to assess the relative 
positions of two vessels. 

Photogrammetric measurements, whose precursor was the 
Frenchman A. Laussedat use the principle of the human eye. 
They use the physical phenomena of perspective and parallax. 
They are based on picture or the individual video frames 
analysis. Photograph (or video frame) is a perspective (middle) 
view of an object on a matrix plane [3], as shown in Fig. 3. 

D f _ Df 
—=— >D= (1) 

On the basis of the known pixel dimensions, camera 
focal length (f) and the real distance between characteristic 
points (D,) it is possible to calculate the distance between 
the recording system and the object according to (1). These 
characteristic points in this vision system are blue and red 
light sources. 

The main aim of the research carried out is to develop an 
algorithm which allows to assess the relative positions of the 
two ships manoeuvring in the vicinity treated as a pair, using 


photogrammetry measurements and basing on structurally 
simple measuring system. The research was carried out on the 
basis of the film realized with the use of construction in which 
LED’s were placed. Construction model was made in 1:2 scale. 
The algorithm has been developed in Matlab environment, 
where also were carried out simulation studies. 


VIEW FROM THE TOP 
Dx 
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Fig. 3. The central perspective projection — representation of the actual 
image on the camera sensor 


MEASUREMENT DEVICES - SYSTEM 
DESCRIPTION 


For carrying out the tests was built measurement system 
consisting of measurement board, rail, camera for image 
acquisition, servo mechanism, laser rangefinder, xPC Target 
(Real-Time System) and notebook. This is shown in the 
Fig. 4. 


Fig. 4. The measurement system 


Measurement board is the place where five light points 
created by LEDs are placed. There are two blue and two 
red points and one green point in the middle, as shown in 
the Figure 5. The board is placed on a trolley on the rail, to 
allow changing the distance between the optical system and 
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the points of light. As a recording system was used Samsung 
S1070 camera with 10MPix CCD image sensor [9]. Its real size 
is 1/2.33”, which in terms of SI units gives the dimensions: 
6.13x4.60mm and the diagonal length is 28mm. The research 
was carried out with a fixed focal length f = 6.3mm of the 
optical system. Video material was recorded with a resolution 
of 640x480[pix]. In order to allow analytical calculations of 
the distance between measurement board and camera image 
sensor, the pixel size p in mm was determined by the formula 
(2). Pixel size estimated from the available technical data of 
the camera [8] is 9.6m. 


(2) 


where: 

— pixel size in [mm], 

length of the longer side of the matrix [mm], 

— length of the shorter side of the matrix [mm], 
number of pixels per longer side of the matrix, 
number of pixels per shorter side of the matrix. 


BLUE LED 


Fig. 5. The measurement board with blue, red and green LED points 


The distance from camera to the measurement board was 
changed by the board with LEDs at small speed movement. 
Camera was in a fixed location. To enable analysis of the 
image located at various angles to the optical axis of the 
camera — recording system was placed on the rotary servo. Its 
rotation angle was controlled by the usage of Real-Time System 
xPC Target. Into the measurement system came also laser 
rangefinder, that due to high accuracy of 103mm was treated 
as a distance pattern. Data from laser rangefinder was recorded 
by real-time system. Program which controlled angle between 
optical axis and measurement board was written in Matlab. 


LED RECOGNITION ALGORITHM FOR VIDEO 
IMAGE 


Photogrammetry is a method for making distance 
calculations on the basis of particularly points which are 
easily distinguishable on the digital image. In the system of 
determining the relative positions of two vessels the board with 
five LED light shown in Figure 5 was used. Light points due 
to the much higher than the average brightness of the image 
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are easily recognizable. The use of light sources of different 
colors allows for easy correlation of points on the image with 
the points in reality. 

The detection of individual points of light is illustrated in 
the diagram shown in Fig. 6. 

The video frame is read, and after that its color space is 
changed from RGB to YCrCb. Through this operation, it is 
possible to analyze in each channel separately, where Y — means 
the total brightness of all pixels in the image, Cr — red difference 
and Cb — blue difference. Color space transformation from RGB 
to YCrCb is given by the formula (3) [1]. 


Y 0.299 0.587 0.114 |R 
Cr|=|-0.169 -0.331 0.5 |G (3) 
Cb 0.5 -0.419 -0.091 B 


The next step is analysis of the image brightness combined 
with the brightest regions selection. Having examined dozens 
of images of LEDs, which were done form several distances 
the distribution of brightness values in Cr and Cb channels for 
blue, green and red LED was verified. On this basis, the mean 
values were established, which became the thresholds in the 
algorithm. By comparing the brightness values around the areas 
selected as the brightest - is assessed whether the bright area 
of the picture is LED. It is also known that the middle points 
of red and blue LEDs (set one above another on the board) 
can not be apart in X direction. If in the picture was found pair 
of desired color diodes — the distance between them is being 
counted. And if not — the next video frame is analyzed. The 
data obtained in the form of LEDs coordinates and the distance 
between pairs of mono-color LEDs, are passed to the algorithm 
for determining the angle and distance. It gives an opportunity 
to determine the distance and angle between the plane in which 
are light diodes and the plane of image sensor. If the pair of 
characteristic points in the image are not found, the algorithm 
leading analytical calculations is not activated. 

This approach helps to eliminate glare, which in the 
case of sunlight are white and in the case of colored light 
(corresponding to the colors of LEDs) — do not fulfill the 
condition of mutual position. 


ALGORITHM DISTANCE AND ANGLE 
BETWEEN LIGHT PLANE AND IMAGE 
SENSOR PLANE FIXING 


Algorithm for distance and angle works on two levels. 
Angle and distance between sensor matrix plane and board 
plane are determined separately on the basis of known real 
distances between LED diodes and corresponding distances 
on the image sensor. In connection with the construction of 
light layout and the occurrence of perspective phenomenon, it 
is known that the length of vertical segments between red and 
blue LEDs are equal only if both planes are parallel. In that case 
the angle between optical axe of camera and light board is equal 
to 90°, which is illustrated in Fig. 3. Analytical calculations 
of distance are not complicated and are carried out according 
to formula (1). There is no need for computing the angle. 
Described case is the simplest one and occurs only in the second 
phase of manoeuver, when two vessels are moving parallel and 
there is no longitudinal shift between them. This part of the 
algorithm enables quick determination of lateral shift, which 
ensures easy and quick identification of the distance between 
the two ships. In the remaining stages of the manoeuver, apart 
the distance is also determined the angle between the optical 
axis of camera and the array of lights. In the algorithm are 
considered four separate cases (see Fig. 7): 
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Fig. 6. LED in image recognition algorithm 
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Fig. 7. a) Mutual LEDs locations on the video image and in reality — case a; b) Mutual LEDs locations on the video image and in reality — case b; 
c) Mutual LEDs locations on the video image and in reality — case c; d) Mutual LEDs locations on the video image and in reality — case d 
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a. LEDs are located on both sides of the matrix center and 
segment between the and the red LEDs is longer than the 
segment between the blue LEDs; 

b. LEDs are located on both sides of the matrix center and 
segment between the and the blue LEDs is longer than the 
segment between the red LEDs; 

c. all points of light are located on the left side of the matrix 
center; 

d. all points of light are located on the right side of the matrix 
center, 

(respectively a, b, c, d). 


The basic dependence used by the algorithm is the 
occurrence of perspective visible in the pictures. The closer to 
the optical system is a segment, that it appears to have a greater 
length. In reality the distance between light points does not 
change. The same rule is used during 3D object from photo 
modeling [11]. Therefore, in order to determine the distance 
between the point where the light is placed (eg. a ship’s bow) 
and the point at which the optical system is set (such as stern 
of the second ship), the intersection points of the image should 
be found. The intersection points are determined by crossings 
of lines contain ing respectively points Al’, B1’ and A1”, B1” 
(Fig. 7). 

Then is determined the length of the segment c= ere 
which corresponds to the distance between these lines 
designated in the middle of the image sensor. On this base 
according to mathematical relationship (1) distance F is 
computed. This is distance between camera lens and virtual 
light table rotation center C. Distance marked as F, in the 
Fig. 8 is the longitudinal distance, which will be passed to the 
control system. 

VIEW FROM THE TOP 


=Al’ 


B1’=B1" 


Fig. 8. Central perspective projection - determination 
of the angle and distance. 


In order to determine the angle between the plane of the 
measurement system, and the plane of the image sensor matrix 
are the basic assumptions in analytical geometry were used [eq. 
(4), (5), (6), (7)]. Angle a is being counted in each algorithm 
step on the basis of equation (8), where the corresponding 
distances are marked in the Fig. 8. 


_Df 
d mP 


(4) 


X 
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G= Fy) (5C tn (6) 
Pa =F- F (7) 
a = arcsin Ta (8) 
C 
where 
D, — the real distance between two red and two blue 
LEDs, 
f — the focal length of the camera lens, 
p  —_ pixel dimension [mm]. 


Distance is being counted in mm and angle a in degrees. 
Algorithm recalculates the angle value into system values 
0° + 180°, where 0° is parallel movement and 180° — means 
that ship is moving parallel but on the other side. The second 
case is almost impossible to obtain in practice. 


RESULTS 


Research was carried out on the basis of the real model of 
the system made on a scale 1:2 (Fig. 4). As the first distance 
assessment algorithm was tested. Calculations, under which was 
designated the distance between the camera and measurement 
table, were made based on the video. It was filmed at the time 
when the light board was moved away parallel to the optical 
axis of the camera along the rail. 

2500 


Laser rangefinder distance 
Distance counted 


5 6 BÈ 


Distance [mm] 


8 


0 20 40 60 60 100 120 
Time [s] 


Fig. 9. Distance count and laser rangefinder measurement. 


During the measurement the distance between the optical 
system and the table was recorded using a laser rangefinder. 
The results are presented in Fig. 9. 

Line described by “laser rangefinder distance” and pointed 
by arrow indicates the reference distance carried out using 
laser rangefinder. The results of computations done in Matlab 
environment using distance evaluation algorithm are also 
presented in Fig. 9 and described “distance counted”. The 
biggest difference between reference and computed distance 
occurs at the end of simulation and reaches 30mm at range 
of 2.13m. It is about 1.55%. At the rest of simulation this 
difference is less than 1%. 
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In the second stage of the study measurements, which 
allowed for evaluation of the correctness of the algorithm 
is the simultaneous assessment the angle and distance, were 
performed. For this purpose the rangefinder combined with 
a camera was placed on the rotary servo. To check if the 
characteristic of servo’s rotation angle is linear, the calibration 
measurement was made. Rotating rangefinder system was 
placed at a fixed distance from the large size vertical calibration 
board. There was performed measurement of angles from 50 
to 130 degrees. Servo rotation angle was set based on the 
basic trigonometric and registered according to laser distance 
measurement, which is described by “laser count” in the 
Fig. 10. 
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Fig. 10. Linearity of the servo rotation characteristics check 


Servo rotation characteristic is close to linear one and for 
the algorithm will be so treated. It could be described by the 
following equation of straight line: f(x) = 0.97x-55.93. 
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Fig. 11. Distance between LEDs estimation 


Array was placed at a fixed distance of 1.52m from the 
optical system. The camera that is treated as a video recorder 
was positioned on the rotary servo motor. Research was carried 
out on the basis of the recorded video in postprocessing mode. 
The control signal angle ofthe servo and data from measuring 
the distance both recorded in real time have been adopted as 
a benchmark. They are described by “servo angle” and “laser 
rangefinder” on the charts presented in Figs 12 and 13. 

The distance between the centers of diodes with the same 
color were calculated by the algorithm and plotted in the Figure 
11 marked with “Ydst”. Similarly is calculated the distance 


between the centers of diodes located in one horizontal line. It is 
marked with “Xdst” in the Figure 11. Based on these distances 
are calculated angles and distances in the subsequent stages of 


the algorithm. 
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Fig. 12. Angle between the plane of light emitting diodes 
and the camera sensor 
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Fig. 13. Distance between the plane of light emitting diodes 
and the camera sensor 


Angle and distance estimations are shown in the graphical 
form in Figs 12 and 13. The biggest difference in between 
reference and counted distance is about 21mm. Which gives 
1.38% of whole distance. These differences are especially 
noticeable at small angles between the planes and when they 
are nearly parallel. The first mentioned case is probably due to 
defects in optics. In the second one occurs difficulty in algorithm 
cases switching. Also the angle difference, which reaches 5° 
at angle of 88°, could be caused by the same process. It gives 
value at the border of acceptability about 5%. At time point of 
60s distance count of 0 mm and angle count of 0 deg occurs. 
There is shown case where 4 characteristic light points were 
not distinguished in the analysed image. In real time system in 
this case samples are being filtered and the distance and angle 
are predicted on the basis of the previous counts. 


CONCLUSIONS 


- Presented algorithm can be successfully used in the control 
system of the relative positions of two vessels moving 
parallel. It is fast and not very complicated. Another positive 
value is the ease of the installation of the cameras and 
matrices of light points on board. Moreover the system is 
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not very sensitive to external factors such as scratches on the 
ship’s sides, rainfall, lack of good lighting and a small mist 
in comparison to system based on laser measurements. 
Errors in position and angle estimation at the level of 1.5% 
of range are acceptable. The results can be improved by 
taking into account the calibration and optical defects of the 
registrant. In the area where errors are the biggest is where 
ships achieve second phase of motion. It may be an area 
for additional laser rangefinders arranged perpendicular to 
the side at the bow and stern. 

Future work includes algorithm running in real time and 
incorporate optical defects such as distortion at short focal 
lengths to improve its performance. 
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